Joachim Schöpfel
email: joachim.schopfel@univ-lille.fr

WHAT ABOUT ODTS? ARE THEY GREY?

Keywords: Theses, dissertations, heritage collections, scientific heritage, grey literature

The term of grey literature is sometimes applied for older material and special collections, especially in the field of digitization projects of scientific heritage. The following paper will analyse this term of "grey scientific heritage" and, based on empirical and conceptual elements, contribute to a better understanding of grey literature. Special attention will be paid on older theses and dissertations (OTDs), as a main part of scientific heritage especially from universities.

Introduction

In systematic literature reviews, meta-analyses and library guidelines, grey literature is often described as unpublished material disseminated outside commercial channels, not peerreviewed and with limited information referencing. The Encyclopaedia of Library and Information Sciences defines grey literature as material "produced on all levels of government, academics, business and industry in print and electronic formats, but which is not controlled by commercial publishers, i.e., where publishing is not the primary activity of the producing body" [START_REF] Schöpfel | Grey literature[END_REF]. The Prague definition of grey literature raises awareness of its documentary nature, intellectual property rights and quality issues [START_REF] Schöpfel | Towards a Prague definition of grey literature[END_REF]. More recently, the 2014 Pisa Declaration on Policy Development for Grey Literature Resources1 mentions inter alia "research and technical reports, briefings and reviews, evaluations, working papers, conference papers, theses and multimedia content" and recognizes them as an "essential resource in scholarly communication, research and policy making, (and) a key source of evidence, argument, innovation and understanding in many disciplines, (and) an important and valuable part of research and information", which generally means that this is current material with new and unpublished results. However, the term grey literature is also applied to older material and special collections, especially in the field of scientific heritage digitization projects. The following paper will analyse the term "grey scientific heritage" and, based on empirical and conceptual elements, contribute towards a better understanding of grey literature. Special attention will be paid to older theses and dissertations (OTDs) as an essential part of scientific heritage, especially from universities. We recently published a debate on whether and why electronic theses and dissertations (ETDs) should still be considered grey literature in the digital age [START_REF] Schöpfel | Are electronic theses and dissertations (still) grey literature in the digital age? A FAIR debate[END_REF].

Here, we will try to assess the grey characteristics of older items from the Gutenberg era. The methodological approach is twofold:

• First, relevant papers on scientific heritage from the last two decades will be reviewed, with special attention paid to the definition of grey literature and to the inclusion of theses and dissertations. • Second, the paper will assess more than one hundred digitization projects from a recent French public digitization program to valorise otherwise hidden scientific heritage collections through digital libraries and open repositories. Do these projects contain grey literature? What kind of grey literature? What role is played by OTDs? This paper will discuss these papers and projects in the light of the usual definitions of grey literature.

Grey literature and heritage -an overview

A search via Google Scholar, OpenGrey and the GreyGuide2 reveals that few papers make explicit links between grey literature and scientific or cultural heritage collections; a textmining study of a large corpus of papers on grey literature confirmed that the bigram "cultural heritage" is rarely used [START_REF] Bartolini | A terminological "journey" in the Grey Literature domain[END_REF]. Most papers on grey literature deal with recent items, not older material.

Scientific vs cultural heritage

The meaning of scientific heritage is also a matter of discussion. The term itself is "diverse, complex, multi-layered (and) difficult to define" (Lourenco & Wilson 2013, p.745). It is part of cultural heritage, i.e. a shared collective legacy, a corpus of material signs handed on by the past in every culture and which constitutes a source of identity and cohesion for communities, everything "we want to keep, share with others and pass on to the next generation" (idem) 3 . Heritage of research is one part of these tangible or intangible assets, everything researchers have produced and what is of interest for future research. In other words, "what the scientific community as a whole perceives as representing its identity, worth being passed on to the next generation of scientists and to the general public as well" (idem, p.746). Yet some ambiguities remain, for instance about the legacy character of research libraries, and if scientific heritage must be a scientific contribution or whether it can also include anything produced by scientists, for instance private materials, diaries, letters etc.

The papers on grey literature that deal with heritage collections are from very different domains, including for example Holocaust literature, urban planning, Polish underground literature, Newton's journal, computer science, Antarctic research, Iceland research publications and the Serbian cultural enlightenment. Some of these papers clearly focus on cultural heritage without any scientific legacy character, such as two studies on historical documents produced by public authorities (de Biagi & Puccinelli 2016) and on Yizkor books [START_REF] Jones | Yizkor Books as Holocaust Grey Literature[END_REF]. Other examples are papers on private collections in the Prado Museum [START_REF] Docampo | Creating a heritage collection: the entry of three private libraries into the Prado Museum Library[END_REF]) and on the Australian Baptist heritage collection [START_REF] Burn | The Australian Baptist Heritage Collection: implications for the management of geographically distributed special collections[END_REF]. These items and holdings may be of interest for scientists but are not produced by scientists.

Other studies include both types of heritage, scientific and not, like a paper on unpublished material (i.e. manuscripts, letters, photographs and sketches) by Sir Julius von Haast, a New Zealand scientist in the 19th Century [START_REF] Nolden | Sir Julius von Haast: exploring an archival documentary heritage collection in the Alexander Turnbull Library[END_REF], an analysis of Polish unpublished and prohibited "underground" literature, including translations of scientific and technical items [START_REF] Nahotko | Some Types of Grey Literature: A Polish Context[END_REF], or the presentation of the Virgin Islands Heritage Collection with a core collection of digitized material with funeral booklets, historical photographs and newspaper articles, but also research reports and occasional papers produced by local research units [START_REF] Marsicek | Virgin Islands Heritage Collection: A Global Resource Sharing Model[END_REF].

Types and definitions of grey scientific heritage

Some papers explicitly address grey literature in terms of scientific heritage, especially scientific and technical reports, working papers, proceedings and surveys (cf. [START_REF] Japzon | Wallops Island Balloon Technology: Can't see the Repository for the Documents[END_REF][START_REF] Stock | LARA-Open access to scientific and technical reports[END_REF][START_REF] Juliusdottir | For better or for worse: Knowledge output 1944-2001 and effects of the Legal Deposit Act no. 20/2002 and e-publishing on access to GL in Iceland[END_REF]. Less frequent grey items in these papers include technical drawings [START_REF] Jackson | Grey Literature and Urban Planning: History and Accessibility[END_REF][START_REF] Biagioni | A Hidden Treasure on Computer Science Pre-History in Pisa: The CSCE Collection[END_REF], newsletters and workshop/training materials (Ramos-Lun & Vogel 2005), handwritten notes (Cirkovic 2015) and materials from conferences that are not readily available [START_REF] Gheen | Digitizing Grey Literature from the Antarctic Bibliography Collection[END_REF]. Two studies mention older theses and dissertations [START_REF] Costello | Louisiana Coastal Wetlands and Louisiana Coastal Grey Literature: Vanishing Treasures[END_REF][START_REF] Biagioni | A Hidden Treasure on Computer Science Pre-History in Pisa: The CSCE Collection[END_REF]. All these items were produced by scientists, are scientific and technical information and written for the scientific community.

Figure 1. Cultural heritage, scientific heritage, grey heritage, grey literature How do these papers define grey literature? Some papers just provide a typology of documents considered as grey. Ramos-Lun & Vogel (2005) cite the GreyNet Luxemburg definition of grey literature and distinguish it from "ephemera", i.e. "materials (physical or electronic) that regardless of appearance, quality or quantity, and that at some point were considered disposable and of little value or no value, through time, had become valuable in such a way that it had broadened their appeal and made them desirable to be collected and preserved by individuals, collectors and information institutions". In fact, they describe "ephemera" as a special kind of grey literature. [START_REF] Rucinski | The Elephant in the Room: Toward a Definition of Grey Legal Literature[END_REF] comments that the ephemeral and variable nature of grey publication types, editions and formats makes them hard to describe and define. The risk of being lost if no investment is made is one of the major characteristics attributed to grey heritage items; they are considered endangered due to their small production quantities and various preservation challenges [START_REF] Jackson | Grey Literature and Urban Planning: History and Accessibility[END_REF]. Another elements commonly used to distinguish grey items from other scientific heritage is the fact that they were not (or not widely) published and may even be confidential, secret or prohibited ("underground"), non-traditional materials, or writings kept away from public review [START_REF] Hofe | Towards a Genealogy of Grey Literature via Newton's Journals[END_REF], Ramos-Lun & Vogel 2005[START_REF] Nahotko | Some Types of Grey Literature: A Polish Context[END_REF]). [START_REF] Juliusdottir | For better or for worse: Knowledge output 1944-2001 and effects of the Legal Deposit Act no. 20/2002 and e-publishing on access to GL in Iceland[END_REF] states that these items are often "neither bibliographically accessible in catalogues open to the public nor available through traditional market publishing distribution channels", while [START_REF] Jackson | Grey Literature and Urban Planning: History and Accessibility[END_REF] adds that they are usually difficult to discover or obtain, often difficult to find in libraries, in online databases and on the web, (and) non-accessible.

The distinction between current items and older material (50 years or more) is not always simple. Some studies are "borderline" and do not make it clear if the main purpose is the discovery of scientific heritage or access to recent resources. For instance, one of the first papers on grey literature remains ambiguous about the wartime German research reports grey literature scientific heritage cultural heritage collected by the US Army and the UK National Lending Service (now British Library Document Supply Centre) from 1945 on [START_REF] Chillag | From Weimar to Maastricht and beyond: half a century with grey literature[END_REF] -was this historical material, or did the US Army gather this material and ship it westward because of its value for cutting edge technology?

Preventing the risk of loss

Why is it important to invest in older grey literature? The papers provide different reasons, e.g. their historical value as "hidden treasures" [START_REF] Stock | LARA-Open access to scientific and technical reports[END_REF][START_REF] Biagioni | A Hidden Treasure on Computer Science Pre-History in Pisa: The CSCE Collection[END_REF], their interest for institutional history and commemoration [START_REF] Anderson | Harnessing NASA Goddard's Grey Literature: The Power of a Repository Framework[END_REF], like the celebration of the 50th anniversary of the first Italian computer [START_REF] Biagioni | A Hidden Treasure on Computer Science Pre-History in Pisa: The CSCE Collection[END_REF]. In other words, for these authors the value is the record of progress, not the information itself.

The value of older materials can be increased by making them more readily available. "The greatest challenge remaining for our library is to make our grey literature and ephemera collection available to our users" (Ramos-Lun & Vogel 2005). A review of older material can include the retrospective enhancement of descriptive and name authority records, thus resulting in "improved documentation of the collections, thereby maximising the discoverability of historical evidence and utilisation of the informational value of a personal and scientific archive" held in library collections [START_REF] Nolden | Sir Julius von Haast: exploring an archival documentary heritage collection in the Alexander Turnbull Library[END_REF]. This means a systematic search and collection of resources from a variety of agencies and organizations [START_REF] Costello | Louisiana Coastal Wetlands and Louisiana Coastal Grey Literature: Vanishing Treasures[END_REF] or digitization and dissemination on the web of collections that were previously fairly inaccessible [START_REF] Gheen | Digitizing Grey Literature from the Antarctic Bibliography Collection[END_REF].

Another related purpose is digital preservation which, alongside dissemination via web servers, prevents grey literature "from moving further toward the black" (Ramos-Lun & Vogel 2005). Here, various authors are in favour of centralized archives, especially of institutional repositories (Jackson 2004, Stock et al;2009[START_REF] Lynch | Updating the Agenda for Academic Libraries and Scholarly Communications[END_REF]) but other solutions may exist.

Standards

Some grey literature, primarily older items, is poorly described in catalogues and databases. Some initiatives insist on the importance of standards to improve the findability and interoperability of heritage collections. Yet, except for a general call for standards, there is no consensus about which kind of standardgeneric or not, etc.should be applied to grey literature. [START_REF] Jackson | Grey Literature and Urban Planning: History and Accessibility[END_REF], for instance, recommends a current and generic cataloguing standard such as AACR2 for the conversion of print resources to accelerate the inclusion of urban planning resources in online databases and catalogues. [START_REF] Kansa | Googling the Grey: Open Data, Web Services, and Semantics[END_REF] promote a field-specific approach to primary data, i.e. "a common and highly abstracted framework for expressing archaeological observations, their descriptive properties and their contextual relationships".

Two metadata standards are generally accepted, the Dublin Core and XML. [START_REF] Anderson | Harnessing NASA Goddard's Grey Literature: The Power of a Repository Framework[END_REF] developed an extension of the DC, called Goddard Core Metadata Element Set for the metadata of resources produced by the NASA Landsat Legacy Project.

In a quite different environment, namely the Cuban Heritage collection at the University of Miami, [START_REF] Baur | The Home Stretch: Developing Automated Solutions for Legacy Container List Data at the Cuban Heritage Collection, University of Miami Libraries[END_REF] apply the Library of Congress Encoded Archival Description format (EAD-XML) for the conversion of legacy search aids in typewritten, MS Word, PDF, HTML and even poorly executed EAD-XML formats. Standardized and valid EAD-XML mark-up, according to the authors, is crucial to provide deep access to and interoperability of their archival description and metadata.

Legal aspects

Some papers address legal issues, but this does not seem to be a major concern for this kind of project. Of course, intellectual property rights must be assessed thoroughly before digitizing heritage collections [START_REF] Blackwell | Hijacking Shared Heritage: Cultural Artifacts and Intellectual Property Rights[END_REF]). Yet the general idea is that older primary source and public domain materials should be freely available, without user agreements and terms of service as a precondition for content access: "(…) our cultural heritage is vulnerable, and risks becoming encompassed within a modern enclosure movement if action is not taken" [START_REF] Clark | Enclosing the public domain: The restriction of public domain books in a digital environment[END_REF]. Licensing and laws should be in favour of projects designed to foster the preservation and accessibility of hidden and rare material which may otherwise be "locked up under copyright" [START_REF] Lynch | Updating the Agenda for Academic Libraries and Scholarly Communications[END_REF].

The case of a national digitization program

As mentioned above, scientific heritage can be described as tangible or intangible assets, everything researchers have produced and what is of interest for future research. The last aspect, in particular, was a central element of a digitization program launched in 2013 and 2014 by the French Ministry of Higher Education and Research to valorise scientific heritage collections4 . The funding criteria were above all the collections' interest for research communities, together with the technical quality, accessibility (open access), interoperability and the added value of the service environment. The document types were not a specific condition of the program, and the participantsmainly academic libraries and documentation centreswere invited to submit proposals including all kinds of library materials such as journals, books, unpublished papers, posters, photographs, maps and so on. Also, the program was not limited to scientific documentsliterary, political or business items could be part of the proposals, as well as press products or personal archives, if their value for research could be established.

123 proposals were submitted and evaluated by an expert group of librarians and scientists. Based on the expert evaluation, the Ministry selected 33 projects with a total grant of €1,296,000 for one or two years. The objective was to valorise otherwise hidden scientific heritage collections through digital libraries and open repositories in an environment of open science and new library services. 38 proposals included grey literature (31%), most of them issued by universities. Nearly all (34) involved consortia with two or more partners; seven consortia incorporated international cooperation with institutions from Italy, Canada, Belgium and Germany. 25) of these 38 proposals are 100% (pure) grey literature projects, with nothing else than grey items like reports, field work, memoirs, minutes or simply "unpublished work". Six projects included dissertations, such as 444 theses in medicine from the University of Lyon (19th and early 20th Centuries), about 360 theses in law from the University of Toulouse (19th Century) and about 1,000 older theses from Bordeaux. The other 13 projects included other documents than grey literature, often non-textual items like maps, drawings and photos. It is difficult to give reliable figures on the overall size of these projects in terms of volumes or pages. A best estimate of these "hidden treasures" is at least 22,500 volumes with 2.9 million pages. The theses and dissertations represented about 8,218 volumes and 1.15 million pages. Even if only a smaller part (15%) of the whole number of projects with grey literature included theses and dissertations, their real share in the grey heritage collections which academic libraries and research laboratories submitted for funding can be assessed as 35%-40%.

Older theses and dissertations

In the French digitization program, older theses and dissertations (ODTs) formed a significant part of the academic collectionsoutside of book and journal holdingsfor which investment, digitization and valorisation through platforms, digital libraries etc. was sought. Yet what exactly does "older" theses and dissertations mean? In the six proposals, one part of the theses and dissertations was written during the 19th and early 20th Centuries, while others are from the last 50 years up to 2005. In other words, one part of these heritage collections is already in the public domain while other items are still protected by intellectual property rights.

A comparison with the figures from a union catalogue (Trove), a search engine (BASE), a portal (DART-Europe) and two repositories (TEL, NTK) reveals some interesting insights, even if these figures should be considered with caution because these tools are not really designed for older items, except the Australian Trove catalogue (table 1). These figures are anything but comprehensive or representative5 , yet illustrate two aspects: on the one hand a non-negligible part of theses in academic catalogues are ODTs, and on the other a small part of them have already been made findable and accessible in the new digital and open science environment via metadata digitization and production. In absolute figures, this small part represents more than 10,000 theses (best estimate).

Table 1 distinguishes between two categories, for two reasons: theses more than 100 years old are probably already in the public domain; theses dating from before 1950 were written by authors who have probably finished their academic (or other professional) careers. Both categories can be considered scientific heritage, and the first category facilitates valorisation through digitization programs.

Age is not the only criteria for scientific heritage. Another, even more important criteria, is the material's quality and value for the research community. Academic theses and dissertations have always been considered as the result of at least three years' original, independent and critical thinking by young scientists, subsequently validated by a commission of senior scholars representing an academic institution and the research community. Depending on the discipline, older theses may still be of interest for today's research. In any case, they represent a unique testimony on the history of science and academic life.

Grey heritage

When does grey literature become scientific heritage? The overview on published studies and papers does not provide a clear indication. Age plays a role, as well as legal status (public domain) and the real usage and interest for research in a given field. In some disciplines like medical and life sciences, scientific documents "expire" and become obsolete more quickly than in others, like mathematics or history. In institutions which switched to ETDs many years ago, media might be another criterion, with all native print theses being considered part of the research heritage; yet this criterion makes no sense where print theses are still accepted.

Grey literature is not easy to define, and the fuzzy term of heritage adds to the confusion.

Nevertheless, the analysis of recent papers and other evidence from digitization programs and discovery tools enable three aspects to be clarified:

1. Grey literature is usually defined as being difficult to get because of limited quantities, non-commercial dissemination etc. Yet the problem of scientific heritage is not acquisition because the items are already in library holdings, research collections etc., but are hidden, invisible and not included in bibliographic records. So the first challenge is not the identification of interesting items somewhere outside of the library or the discovery of external information sources and channelsthe real challenge is the discovery of valuable material inside the library or research institutions, and the decision to invest into "digging up the treasure" and make it visible which, in the open science environment, means digitization and online publishing. 2. Grey literature is above all a challenge for acquisition and collection building. Yet the secondary challenge with scientific heritage is not collection building but how to improve the findability and accessibility of formerly hidden items in a new service environment on the internet. Findability means that all items must be described with rich metadata and, if possible, linked to the semantic web, applying the usual standard formats. Accessibility means that these items should be freely available on the web via digital libraries, open repositories or similar platforms, without unnecessary restrictions and without artificial enclosure, and accessible via standard protocols such as OAI-PMH. 3. Grey literature needs curation and conservation efforts. The third challenge for scientific heritage is similar, even if part of the hidden treasures are already well preserved -perhaps too well. In the era of open science, the preservation of scientific heritage means digital preservation in a secured environment but, at the same time, reusability with new digital tools and services, including content mining and linking to data (cf. France 2018).

Scientific heritage like ODTs, as valuable and useful textual sources of information and insofar as it needs discovery, curation and preservation, can thus be considered grey literature. The intermediary role and importance of academic and research libraries remains intact, and their skills, investment and engagement are needed to discover these hidden treasures and make them findable, accessible and reusable. However, the direction of intervention and perimeter of action are shifted from outside to inside, as these treasures are already in the holdings, waiting for discovery on bookshelves or in microfilm containers like Sleeping Beauties.

Conclusion

Grey literature is not easy to define, and the terms cultural and scientific heritage add even more fuzziness to the concept. Our findings can be summarized in three sentences:

• ETDs are scientific grey literature.

• OTDs are grey scientific heritage.

• Some ODTs are less grey than others. 10

There is growing awareness of the economic and social benefits of cultural heritage 6 . Scientific heritage is part of this cultural heritageit is produced by scientists and is of interest for the research community. For over 20 years, theses and dissertations, "the most useful kinds of invisible scholarship and the most invisible kinds of useful scholarship" [START_REF] Suber | Open access[END_REF]) are part of ETD programs, and an increasing number of universities and countries such as Australia, New Zealand and France, have gone 100% digital for theses. It is time to take care of the older theses and dissertations present in many academic and research libraries. The problem with ODTs is not dissemination and acquisition but findability, accessibility and reusability. They must be digitized, described and republished in a modern service environment. We recently asked the question of whether ETDs are still grey literature or not, and concluded that "if by 2020 ETDs should be completely integrated in the emerging open science infrastructures, as open as possible (and just as closed as necessary), easily retrievable and accessible, and largely reusable by content mining tools, greyness would no longer be a problem" [START_REF] Schöpfel | Are electronic theses and dissertations (still) grey literature in the digital age? A FAIR debate[END_REF].

The same conclusion applies to ODTs. Obviously, as they are part of scientific output, intellectual work and valuable for scientific and historical research, they are grey literature because of their limited availability, lack of description and risk of loss. Just like all grey literature, ODTs need care and curation by information professionals, especially in the academic and research libraries in the frontline, to increase their findability, accessibility and potential reuse (figure 3). The challenge of grey scientific heritage is conservation and dissemination through open access for the scientific community. However, in the digital age and the emerging open scientific age, it must be insisted that readers include not only scientists but also interested citizens, and that they are not only humans but also machines, with a corresponding impact on decisions on how to publish and how to describe the documents.

6 See the 2015 report of the European Commission, cf. footnote 3. As the number of older theses and dissertations is by definition limited, we may be hopeful that one day all or most of these documents will be searchable and available on open repositories or via academic portals as an essential contribution to global scientific heritage.

Figure 2 .

 2 Figure 2. The share of grey literature projects in the French digitization program (N=123)Two thirds (25) of these 38 proposals are 100% (pure) grey literature projects, with nothing else than grey items like reports, field work, memoirs, minutes or simply "unpublished work". Six projects included dissertations, such as 444 theses in medicine from the University of Lyon (19th and early 20th Centuries), about 360 theses in law from the University of Toulouse (19th Century) and about 1,000 older theses from Bordeaux. The other 13 projects included other documents than grey literature, often non-textual items like maps, drawings and photos. It is difficult to give reliable figures on the overall size of these projects in terms of volumes or pages. A best estimate of these "hidden treasures" is at least 22,500 volumes with 2.9 million pages. The theses and dissertations represented about 8,218 volumes and 1.15 million pages. Even if only a smaller part (15%) of the whole number of projects with grey literature included theses and dissertations, their real share in the grey heritage collections which academic libraries and research laboratories submitted for funding can be assessed as 35%-40%.

Figure 3 .

 3 Figure 3. Challenges of grey scientific heritage

Table 1 .

 1 Share of older theses and dissertations (in %)

	Online service	Up to 1950	Up to 1900
	Trove	6.4%	0.5%
	BASE	1.5%	0.7%
	DART-Europe	1.0%	0.4%
	TEL	0.5%	0.0%
	NTK GL Repository	0.2%	0.0%

http://greyguide.isti.cnr.it/pisa-declaration/ 11th Conference on Grey Literature and Repositories: proceedings [online]. Prague: National Library of Technology,

 [cit. 2018-12-15]. Available from: http://nrgl.techlib.cz/conference/conference-proceedings/.ISSN 2336-5021.

GreyGuide Repository and Guide to Good Practices and Resources in Grey Literature http://greyguide.isti.cnr.it/

Cf. also https://en.unesco.org/ By the way, the report of the Horizon 2020 expert group on cultural heritage commissioned by the European Commission and published in 2015 (Getting cultural heritage to work for Europe) simply doesn't define the term.11th Conference on Grey Literature and Repositories: proceedings [online]. Prague: National Library ofTechnology, 2018Technology, [cit. 2018-12-15]-12-15]. Available from: http://nrgl.techlib.cz/conference/conference-proceedings/.ISSN 2336-5021.

11th Conference on Grey Literature and Repositories: proceedings [online]. Prague: National Library ofTechnology, 2018Technology, [cit. 2018-12-15]-12-15]. Available from: http://nrgl.techlib.cz/conference/conference-proceedings/.ISSN 2336-5021.

The digitization program was part (segment

5) of the Bibliothèque Scientifique Numérique (BSN) or Digital Scientific Library framework(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). The corresponding author was a member of the steering committee in charge of project selection. The call, criteria and results of the program are available at http://www.bibliothequescientifiquenumerique.fr/bsn-5-numerisation/

For instance, we could have added the French academic union catalogue SUDOC or the catalogue of the National Library of the Czech Republic; but would the results have been very different?11th Conference on Grey Literature and Repositories: proceedings [online]. Prague: National Library ofTechnology, 2018Technology, [cit. 2018-12-15]-12-15]. Available from: http://nrgl.techlib.cz/conference/conference-proceedings/.ISSN 2336-5021.


J J Da Costa

F Chainet

B Celse

M Lacoue-Nègre

C Ruckebusch

D Espinat

Comparing kriging, spline and MLR in product properties modelization

Keywords: Kriging, Multilinear Regression, Splines, Interpolation, Hydrocracking, Diesel, Cloud point 27 28

In the most part of chemometrics application concerning the prediction of physico-chemical properties, 12 regression models are preferred as they are easy to implement and their posterior analysis provide simple 13 interpretations. Interpolation methods are most currently used in geostatistics applications where the dimension 14 of study area is generally limited. In this work, we proposed to develop kriging or splines models for predicting 15 petroleum products properties. Kriging and splines have different foundations as the former is based on 16 stochastic assumptions and the latter is built on deterministic approach. A well-illustrated comparison of both 17 methods was carried out through three suitable examples to point their similarities and their divergences. The 18 advantages of using kriging or / and splines instead of classical regression models were also discussed. Results

19

pointed out the flexibility of interpolation methods as they provide good accuracy for linear and nonlinear cases.

20

They also confirmed previous studies which pointed out equivalence between kriging and spline models 21 performances in some situations. However, kriging approach has more valuable aspect among other interpolation 22 methods since it provides a measure of prediction uncertainties. Kriging modeling were finally compared to 23 multilinear regression for the prediction of diesel cloud point ranging from -39°C to -12°C. Models 24 performances pointed out that kriging enables to improve both accuracy and robustness.

Material and methods

20

All the datasets used to test predictive methods in the present work are described here. Main theoretical aspects 21 of the used methods are then recalled.

22

Parametric regression models 23

In this section, ݖ denotes a variable to model and ݔ a set of explanatory variables. It is assumed that a set of 24 measured values ݕ( ଵ , … , ݕ ) at ݊ sampled data ݔ( ଵ , … , ݔ ) was observed. In addition, all measured values are 25 considered to be submitted to an error due to experimental process such that:

26

= ࢠ + ࣕ , = , … , (Eq. 1)

The main idea of parametric regression models is to approximate ݖ by a known mathematical function ݂ such 1 that [START_REF] Hastie | The elements of statistical learning: Data mining, inference, and prediction[END_REF]):

2

= ሺ , ሻ + (Eq. 2)

3
where ߠ denotes the set of parameters related to the analytic expression of ݂ and ߝ is the modelling error at the 4 ݅ ௧ sampled data. In what follows, ߠ is assumed to be a -dimensional vector (݊ >  > 1).

5

Multilinear regression 6

In case of multilinear regression (MLR), ݂ is assumed to have the following form [START_REF] Azaïs | Le modèle linéaire par l'exemple: Régression, analyse de la variance et plans d'expérience illustrés avec R et SAS[END_REF]:

7 ሺ, ሻ = + (Eq. 3)
8

where ݔ ் denotes the transpose vector of .ݔ Let ܺ denotes the design matrix i.e. the matrix whose the ݅ ௧ row 9 contains the coordinates of the ݅ ௧ sampled point, and ݕ = ሺݕ ଵ , … , ݕ ሻ ் . The goal of MLR is to estimate "as well 10 as possible" the set of parameters ߠ. The most current approach consists on introducing Gaussian assumptions 11 for modeling errors. The ߝ are then generally assumed as independent identically normally distributed random 12 variables with mean zero value and of unknown variance ߪ ଶ . Under these assumptions, the least of squares (LS)

13

estimators which consists on minimizing the sum of squared deviation has optimal characteristics according to 14 inferential statistics [START_REF] Azaïs | Le modèle linéaire par l'exemple: Régression, analyse de la variance et plans d'expérience illustrés avec R et SAS[END_REF][START_REF] Hastie | The elements of statistical learning: Data mining, inference, and prediction[END_REF]). Furthermore, it has been shown that for MLR, the LS-estimator ߠ ௌ can be expressed as a linear combination of measured values, [START_REF] Hastie | The elements of statistical learning: Data mining, inference, and prediction[END_REF]):

16

= ሺ ሻ ି (Eq. 4)

17

This result may be extend to generalized multilinear regression models [START_REF] Azaïs | Le modèle linéaire par l'exemple: Régression, analyse de la variance et plans d'expérience illustrés avec R et SAS[END_REF]Dror and 18 Steinberg 2008). Let ݃ denotes a mathematical differentiable function such that:

19

ሺ, ሻ = + ሺሻ (Eq. 5)

20

then, the LS-estimator has the following generalized expression [START_REF] Hastie | The elements of statistical learning: Data mining, inference, and prediction[END_REF]):

21

= ሺሺሻ ሺሻሻ ି ሺሻ (Eq. 6)

22

where ݃ሺܺሻ denotes the matrix of ݔ by the transformation ݃.

23

Although the LS-estimator calculation results in a relatively simple matrix multiplication, it is rigorously 24 based on deterministic approach which does not include the measurement error. Thus, other estimators related to 25 the concept of regularization were developed in statistical learning (Ridge, Lasso, etc.) (Bolasso: Model The main advantages of using regression models is that the relation between input variables and the 1 modeled response are obvious to interpret. However, find the adapted form of the regression function ݂ may be 2 very complex, especially in nonlinear case.

3

Interpolation methods

4

All the definitions made in introduction of the previous section are still effective in what follows. The basic idea

5

of interpolation is to estimate the value of a variable at a given point within the study area by computing a 6 weighted average of the observed value [START_REF] Arnaud | Estimation et interpolation spatiale: Méthodes déterministes et méthodes géostatiques[END_REF][START_REF] Cressie | Statistics for spatial data[END_REF]). Interpolation methods can be 7 divided into two groups:

8

• deterministic method where ݖሺݔሻ is considered as a "regionalized variable" [START_REF] Wackernagel | Multivariate geostatistics: An introduction with applications[END_REF],

9

• probabilistic method where ݖሺݔሻ is assumed to be a realization of random variable [START_REF] Goovaerts | Geostatics for natural Resources evaluation[END_REF].

10

Both types of interpolation methods are illustrated below through spline and kriging approaches.

11

Spline interpolation 12

The aim of splines interpolation is to draw contour lines "as smooth as possible", that is a map which

13

looks like what a draftsman would obtain manually [START_REF] Dubrule | Comparing splines and kriging[END_REF]. The splines consist of polynomials of degree 14 ݉ being local. The polynomials describe pieces of a line or surface. For degree ݉ = 1, 2 or 3, a spline is called linear, quadratic or cubic respectively [START_REF] Ahlberg | The Theory of Splines and Their Applications[END_REF]. In order to simplify the theoretical aspect of spline 16 the study area is assumed to be a subset of the two-dimensional space. Let ሺ,ݑ ݒሻ denotes the coordinates of .ݔ which is solution of the following minimization problem:

19 ܖܑܕ ࣘ ሺࣘሻ = ∬ ቈቀ ࣔ ࣘ ࣔ ቁ + ቀ ࣔ ࣘ ࣔ ቁ + ቀ ࣔ ࣘ ࣔࣔ ቁ ℝ (Eq. 7)
20 among all the functions ߶ honoring the sampled data points. ܬ ሺ߶ሻ can be interpreted, in first approximation as 21 the bending energy of a thin plate represented by ߶ [START_REF] Dubrule | Comparing splines and kriging[END_REF]. Thus, the function that minimizes ܬ ሺ߶ሻ 22 takes the shape of a thin plate which would be forced to pass through sampled data points. However, this 23 approach does not consider the measurement error related to the observed data. In practice, a most useful approach called "smoothing spline" is preferred. It consists on adding a term of regularization to the objective 25 function. The function to minimize become [START_REF] Dubrule | Comparing splines and kriging[END_REF][START_REF] Wahba | Some New mathematical -Methods for Variational Objective Analysis Using Splines and Cross Validation[END_REF]:

26 ∑ ୀ ሾࣘሺ ሻ -ሿ + ሺࣘሻ (Eq. 8)
ߣ is called the "smoothing parameter" and the weights ݓ ଶ are generally chosen inversely proportional to the 1 error variance. Consequently, the added term force the solution to pass "not too far" from the measured values 2 [START_REF] Dubrule | Comparing splines and kriging[END_REF]. [START_REF] Duchon | Spline associated to n observations of a random function[END_REF] showed that interpolating and smoothing splines have the same 3 expression. For the two-dimensional case:

4

ࢠ ො ሺ , ሻ = + + + ∑ ୀ ܗܔ (Eq. 9)

5

where ݎ represents the Euclidian distance between ݔ and ݔ , i.e.:

6 = ሺ -ሻ + ሺ -ሻ (Eq. 10)
7

However, their coefficients ܽ , ܽ ଵ , ܽ ଶ and the ܾ are obviously different. In case of smoothing spline they are 8 obtained by solving the following system: (Eq. 11) 10

9 ൞ ∑ ୀ = ∑ ୀ = ∑ ୀ = ࢠ ො ሺ , ሻ + = ሺ , ሻ, ∈ ሼ, … . , ሽ,
These results may be extended to -dimensional case. The generalized expression of ܬ ሺ߶ሻ is given below 11 [START_REF] Wahba | Some New mathematical -Methods for Variational Objective Analysis Using Splines and Cross Validation[END_REF]:

12 ሺࣘሻ = ∑ ! ࢻ !ࢻ !…ࢻ !  ⋯  ቆ ࣔ ࣘ ࣔ ࢻ …ࣔ ࢻ ቇ … ℝ ࢻ ା⋯ାࢻ ୀ (Eq. 12)
13

The solution of this generalized minimization problem is discussed in [START_REF] Wahba | Some New mathematical -Methods for Variational Objective Analysis Using Splines and Cross Validation[END_REF].

14

Developing smoothing spline model require to find a suitable value for both smoothing parameter and 15 polynomial order. In practice, a statistical method called generalized cross-validation is usually computed in the 16 learning step [START_REF] Wahba | Some New mathematical -Methods for Variational Objective Analysis Using Splines and Cross Validation[END_REF]. 17

Kriging modeling 18

The goal of kriging is to obtain the best linear unbiased estimate (B.L.U.E.) according to the observed data. Let

19

̂ݖௌ ሺݔ ሻ denotes the kriging estimate of the response ݖ at a given point ݔ . Thus, it must satisfy the following 20 assumptions:

21 ቐ ሺሻ ࢠ ො ሺ ሻ = ∑ ሺ ሻ ୀ ሺሻ ሾࢠ ො ሺ ሻሿ = ࢠሺ ሻ ሺሻ ࢠ ො ሺ ሻ = ܚ܉ ܖܑܕ ࢠ ො ሺ ሻ ሾࢠ ො ሺ ሻ -ࢠሺ ሻሿ (Eq. 13)

22

where ̂ݖሺݔ ሻ denotes any linear estimator of ݔሺݖ ሻ. The definition above calls for a predefined probabilistic the set of values ሼݖሺݔሻ, ݔ ∈ Αሽ describes a multiGaussian stochastic process [START_REF] Stein | Gaussian approximations to conditional distributions for multi-Gaussian processes[END_REF]. According to the 25 process characteristics, two types of kriging may be distinguished. The method is called "simple or ordinary" kriging when the stochastic process is assumed to be stationary of order two, i.e. the expected value is invariant 1 within the study area and the covariance between two measured values only depend on the deviation vector of 2 their related locations [START_REF] Goovaerts | Geostatics for natural Resources evaluation[END_REF]. Otherwise, the method is called "universal" and the expected value has 3 the analytical expression of a polynomials [START_REF] Dubrule | Comparing splines and kriging[END_REF]. 4

Simple kriging predictor 5

The multiGaussian assumptions related to simple kriging results in the two following 6 constraints:

7 ቊ ሾሺሻሿ = ൫ࢠሺ ሻ, ࢠሺ ൯ሻ = ሺ -ሻ (Eq. 14)
8

where ܧሾݕሺݔሻሿ is the expected value of ݕሺݔሻ. The conditions of optimality verified by the kriging predictor lead 9

to solve the following linear system [START_REF] Goovaerts | Geostatics for natural Resources evaluation[END_REF]:

10 ൛∑ ሾ൫ -൯ + , ࣕ ൫ ൯ሿ ୀ = ൫ -൯, ∈ ሼ, … , ሽ (

15

This result implies that kriging estimator can be written as a linear combination of the ሼܥሺݔ -ݔ ሻ, ݅ = 1, … , ݊ሽ,

16

i.e. [START_REF] Matheron | The theory of regionalized variables and its applications[END_REF]:

17 ࢠ ො ሺ ሻ = ∑ ሺ -ሻ ୀ (Eq. 17)
18

where the coefficients ܾ are solution of the following system:

19 ቊ ∑ ୀ = ࢠ ො ൫ ൯ + ࣕ ሺ ሻ = ൫ ൯, ∈ ሼ, … , ሽ (Eq. 18)

20

The expression above is analogous to Eq.11 which is related to splines interpolation. Thus, the kriging predictor 21 may be interpreted as an exact interpolator.

22

Universal kriging predictor 23

A more generalized kind of kriging so called "universal kriging (UK)" or "kriging regression" was also used in 24 the literature [START_REF] Dubrule | Comparing splines and kriging[END_REF][START_REF] Goovaerts | Geostatics for natural Resources evaluation[END_REF][START_REF] Matheron | Splines and Kriging: Their Formal Equivalence: Fontainebleau[END_REF]. It consists on introducing a polynomial trend 25 depending on the available observations, such that the may be expressed as follows:

26 ࢠሺሻ = ∑ ࢼ ሺሻ ୀ + ࢻሺሻ (Eq. 19)
where ݉ denotes the degree of the polynomial trend and ߙሺݔሻ is the approximation error. The analytical 1 expression of components ݂ ሺݔሻ are given in Appendix C. It has been shown that the UK weights predictor are 2 solution of the following system:

3 ቊ ∑ ሺ ሻ ୀ = , = , … , ࢠ ො ൫ ൯ + ࣕ ሺ ሻ = ൫ ൯, ∈ ሼ, … , ሽ (Eq. 20)
4 5

All the system related to the calculation of kriging weights predictor depends obviously on the chosen 6 covariance model .ܥ Some examples of covariance function are given in Table 1. 7

ߪ ଶ ݁ ି ∑ ఏ ೕ ห௦ ೕ ି௪ ೕ ห ೕసభ Gaussian ߪ ଶ ݁ ି ∑ ఏ ೕ ሺ௦ ೕ ି௪ ೕ ሻ మ ೕసభ Linear ݉ܽݔ൛0,1 -ߠ หݏ -ݓ หൟ Spherical 1-1.5ߦ + 0.5ߦ ଷ , ߦ = min ൛1, ߠ หݏ -ݓ หൟ 9
In the simplest case, the measurement errors variance is assumed to be homogeneous around the study area. In 10 case of a heterogeneous errors variance the method is then called "cokriging" [START_REF] Dubrule | Comparing splines and kriging[END_REF][START_REF] Goovaerts | Geostatics for natural Resources evaluation[END_REF].

11

Kriging also requires previous analysis to choose the most adapted covariance structure. Geostatisticians use a

12

graphical recognition tool called "semi-variogram" [START_REF] Goovaerts | Geostatics for natural Resources evaluation[END_REF][START_REF] Isaaks | Applied geostatistics[END_REF]). However,

13

there is no method that enables to compute this step automatically. In this work we only used a Gaussian model

14

(see Error! Reference source not found.) for all developed kriging models. Thus, the learning step consists on 15 estimating intrinsic parameters of Gaussian covariance (including a possible nugget effect, see Appendix E). 16

This step is usually computed by maximizing the likelihood function (see Appendix C) (Isaaks and Srivastava

18

In what follows, only simple kriging was considered.

19

Prediction uncertainty 20

Most of the predictive methods are deterministic as they have no assessment of errors with predicted values. The 21 traditional approach for modeling prediction uncertainty at an unsampled location ݔ consists on computing a minimum error variance ݕ ොሺݔ ሻ of the unknown value ݔሺݕ ሻ and the associated variance ߪ ଶ ሺݔ ሻ = ݕሾݎܸܽ ොሺݔ ሻ -1 ݔሺݕ ሻሿ. The estimate and error variance are then typically combined to derive a Gaussian-type confidence 2 interval centered on the estimate value [START_REF] Goovaerts | Geostatics for natural Resources evaluation[END_REF]. For example, the 95% confidence interval is taken as:

3 4 ሼሺ ሻ ∈ ሾ ෝሺ ሻ -ሺ ሻ, ෝሺ ሻ + ሺ ሻሿ ሽ = . ૢ
(Eq. 21)

5

For example, it has been shown that an estimate of the minimum error variance related to MLR model is 6 expressed as follows (Azaïs and Bardet 2012):

7 ෝ ሺ ሻ = ෝ ሺ + ᇱ ሺ ሻ ି ሻ (Eq. 22)
8

where ߧ ො ଶ is the commonly used estimate of the modeling errors variance which is defined as (Azaïs and Bardet 9 2012):

10 ෝ = ି ∑ ሺ -ෝ ሻ ୀ (Eq. 23)
11

A more rigorous approach consists on assessing first the uncertainty about the unknown, then deduce an estimate 12 optimal in some appropriate sense [START_REF] Journel | When do you need a trend model in Kriging ?[END_REF][START_REF] Srivastava | Maximum variance of profitability[END_REF]. Let ݔሺݑ ሻ denotes the random where ܽ and ܾ are scalars. Note that these probability intervals are independent of any particular estimate ݕ ොሺݔ ሻ

18

of the unknown value ݔሺݕ ሻ, but only depend on the available information ሺܫሻ and ݔ . Each conditional

19

probability distribution function provides a measure of local uncertainty in that it relates to a specific location ݔ .

20

Under the multiGaussian assumptions related to simple kriging (paragraph 2.2.2.1), it has been shown that the 21 mean and the variance of the conditional distribution function at ݔ are equal to the simple kriging estimate

22

̂ݖௌ ሺݔ ሻ and simple kriging variance obtained from the available sampled data [START_REF] Journel | Mining geostatistics[END_REF].

23

2.4 Data sets 24

Simulated data 25

Dataset 1 and dataset 2 are represented in Fig. 1a and Fig. 1b respectively. Each dataset contains 50 26 points uniformly distributed around ሾ-3; 3ሿ ଶ . For the dataset 1, the function to model is affine. For the dataset 2, 27 the function to model is rational fraction. Their analytical expression is given in Table 2. 28 1

Where ߝ were chosen as a normally distribute random variable with mean zero and a fixed variance ߪ ଶ . Details 2 are given in Table 2.

3

In the third situation, we considered a function ݂ which may be interpreted as the sum of two potential 4 functions centered in ܲ ଵ = ሺ-0.5; -0.5ሻ and ܲ ଶ = ሺ0.5; 0.5ሻ respectively. The sampled points were generated 5 such that they describe two well-identified clusters. Then, the coordinates of points were assumed to follow a 6 two-dimensional Gaussian distribution (see Appendix B). Half of coordinates were distributed such that:

7 ሺ , ሻ ~ घ ቆ൫ . . ൯; ቀ . . ቁቇ (Eq. 26)
8

The second part were distributed such that: 8

Models performances were then evaluated using the following process: 9

1. the study area was discretized as a 100x100 grid points regularly spaced to use them as validation set 10 2. both theoretical and predicted values were calculated at each grid point.

11

3. the root mean square error )ܧܵܯܴ( of prediction and the mean relative deviation )ܦܴܯ( were 12 calculated on validation set as detailed in Table 5.

13

For real data, models were compared in two ways:

15

All models performances are discussed in the following paragraphs.

Models comparison on simulated database 1

In each case of simulated data, classical regression, kriging and spline models were developed. The intrinsic 2 parameters estimates related to spline and kriging model (introduced in paragraph Error! Reference source not 3 found.) are specified in Table 6. The models performances are discussed in three ways:

4

• First by overlapping the response surface that the model provided with the original graph of the 5 function to fit.

6

• Secondly by plotting a 2D-graphics of the prediction errors distribution within the study area.

7

• Thirdly by comparing their statistics given in Table 7.

8 9

4

Results for kriging and splines models are illustrated in Fig. 5c and Fig. 5d respectively where the 5 corresponding prediction errors are plotted. The proximity between kriging predictions and real values is also 6 effective. However, some differences may be noted relative to MLR model: Firstly, it may be observed in Fig. 5c 7 that the points where prediction error is close to 0 describe a nonlinear curve which is the intersection line 8 between kriged surface and original map. Then, kriged surface is curved and the quality of the fit is less effective 9 for kriging than for MLR. Indeed, red regions that refers to absolute deviation around 1.5 may be observed in 10 Fig. 5c whereas the maximum value of absolute deviation for MLR model is less than 1 (see Fig. 5b).

11

The curvy characteristic of the response surface may be extended to the thin plate spline (Fig. 5d).

12 However, the prediction error distribution that spline model provides is closer to MLR than kriging. Globally, all 13 obtained models provided good performances since blue regions are predominant in absolute deviation plots 14 (Fig. 5b, c andd). Models statistics given in Table 7 are in accordance with these observations above. Indeed, a 15 MRD value of 0.07 was obtained for MLR against 0.10 and 0.13 for spline and kriging respectively.

Data set 2 (rational function) 7

For simulated data set 2, the function to fit is a rational fraction. Therefore, MLR is obviously low 8 performing here since it is clearly impossible to correctly fit the sampled data by a plane surface. Moreover, find to spline model (Fig. 6b). This relative equivalence between both models is supported by statistics since RMSE 2 value of 1.69 was obtained for spline model against 1.88 for kriging (see Table 7). Note also that the high values For this third example, the same process was used to compare models. As in previous case, we 11 obviously considered that MLR cannot be appropriate to the simulated data observed in Fig. 2. Then, only 12 kriging and spline models performances were illustrated in Fig. 7a and Fig. 7b respectively. As in the previous

Uncertainties of prediction 8

Uncertainties of predicted values were computed by estimating their 95% level confidence interval amplitudes.

9

Results are illustrated for datasets 2 and 3 in Fig. 8

2

Interval amplitude that stochastic approach provided are illustrated in Fig. 8 for dataset 2 and 3. 3 Amplitudes level are given through a color scale (blue regions correspond to lowest interval amplitude and red 4 regions indicate highest amplitudes). It may be observed in Fig. 8a and Fig. 8b Energy & Fuels However, previous examples showed that this polynomial structure is not always appropriate, especially in case 1 of potential function. An advantage in using Gaussian type is that the predicted value always aims at an expected 2 value.

3

Globally, the results showed the adaptability of both kriging and spline models to linear or nonlinear situations 4 whereas classical regression require to find the suitable model structure according to the observed data. Except 5 for classical functions (linear, affine, second degree polynomial, etc.), this may be very complex.

6

Kriging and spline provided very close performances in the tested situations. In addition, the stochastic approach dimensional study remains always challenging. The real difficulty consists on recognizing the adapted 1 covariance structure. On this point, Gaussian structure appears as a good default alternative.

26

  Consistent Lasso Estimation Through the Bootstrap; Hastie et al. 2001; Mohri, Rostamizadeh and Talwalkar 27 2012).

  13 variable modeling the uncertainty about ݔሺݕ ሻ. The distribution function ݔሺܨ , ݖሻ = ݔሺݑܾሼݎܲ ሻ ≤ |ݖሺܫሻሽ made 14 conditional to the available information ሺܫሻ fully models that uncertainty in the sense that probability intervals 15 can be derived, such as: 16 ሼሺ ሻ ∈ ሾ; ሿ|ሺሻሽ = ሺ , ሻ -ሺ , ሻ (Eq. 24)

Fig. 1 5 Fig. 2 :

 152 Fig. 1 : a) 2D plots of dataset 1 in the corresponding study area; b) 2D plots of dataset 2 in the

8 2 .Page

 2 2D points were then randomly generated around the study area as sampled data 9 3. A predefined bivariate function ݂ was applied to each sampled data and a "white noise" (i.e. normally 10 distributed random variable with mean zero and a fixed variance) was added to the result to obtain the 11 corresponding observed responses.12 Thus, each data set consists on ݊ triplet ሺ,ݔ ,ݕ ݖሻ such as 13

Page

  Fig. 3 : Block flow diagram for data acquisition

7Page

  samples) sets. Training samples were selected using a space filling algorithm(Santiago, Sergent 2012) in order to ensure that most predictive situations may be rigorously considered as interpolation 9 cases. The projection of training samples in Conversion-T 95 space is plotted in Fig.4. The color scale is related 10 to measured CP (using NF ISO 23015 (NF EN ISO 23015 1994)). Samples were numbered and labeled from 1 to 11 40. A gradient of CP may be globally observed since CP tends to increase with the T 95 value. In contrast, the 12 higher the conversion rate, the lower the CP is. These observations are in accordance with previous study which 13 pointed that diesel CP and T 95 values are both strongly dependent on its heavy n-alkanes contents. Furthermore, 1 hydrocracking enables to convert molecules combining cracking and isomerization reactions. Then, n-alkanes 2 are mostly converted into branched alkanes when the conversion rate is high. The third variable used to model 3 CP is related to Feedstock which may also significantly affect the diesel quality.

7

  Fig. 4 : Projection of training samples in conversion rate-T95 space

  leave-one-out (LOO) cross-validation was first applied to training data and a ܧܵܯܴ of cross 15 validation and the percentage of predicted points with a precision less than confidence interval (CI) of 16 the standard measurement method were then deduced (see formulae in Table5) (Wahba and models were developed for each simulated database introduced in paragraph 2.4. Their 14 performances were evaluated over the study area and were compared with those of classical regression model.

  13 obviously provides very good performances as all the assumptions of gaussian model are verified. Results are 14 illustrated in Fig. 5. Consider first the results for MLR model. Fig. 5a illustrates the overlapping of the MLR 15 response surface (plane colored in blue) to the original graph (plane colored in black); green points represent the 16 training data. The corresponding prediction errors are plotted in Fig. 5b; the study area is delimited by the two-17 dimensional space and the absolute deviation related to each point is represented by a color scale; white circles 18 localize training points It may be observed in Fig. 5a that MLR response surface and original graph are 19 obviously very close. Note also that blue regions which refers to low prediction errors are predominant in Fig.

  20

  5b. Both above observations show the good quality of the fit, in accordance with the good properties of the LS 21 estimator in case of Gaussian linear model. However, MLR response surface and original graph are not entirely equivalent. Indeed, the visual contrast between blue and black color that appears in Fig.5acombined with the 1 symmetric distribution of prediction errors observed in Fig.5bsuggest a slight tilt of MLR response surface 2 relative to the original plane. Note also that the symmetrical aspect of the prediction errors distribution does not 3 depend on the training set configuration.

Fig. 5 :

 5 Fig. 5 : a)-Overlapping of MLR response surface (blue) and original map (black) for affine function (data

9

  the adapted analytical expression to compute classical regression (linear or nonlinear) appears unlikely. Thus, 10 only kriging and spline response surface are discussed below. Prediction errors distribution is represented in Fig.116a for kriging model. The model appears to provide good accuracy in the most part of the study area. The dark 12 blue regions that refer to the lowest prediction errors are localized around the training samples locations (white 13 circles). A low quality of fit is particularly observed in a well-localized region that correspond to both low 14 samples and a high variability of the function to fit. These observations may be extended 1

3 5 Fig

 5 Fig. 6 a)-Overlapping of kriged response surface (blue) and original map (black) for rational function

  Fig. 7 : a)-Overlapping of spline response surface (blue) and original map (black) for sum of potential

  for stochastic approach. Classical intervals have obviously 10 Page 20 of 30 ACS Paragon Plus Environment Energy & Fuels ellipsoidal form due to their analytical expression. The ellipsoid size increases regularly with the distance from 1 global barycenter of training data.

Fig. 8 :

 8 Fig. 8 : Amplitude of 95% level confidence interval within the study area; a) for dataset 2; b) for dataset 3

7Fig. 9

 9 Fig. 9 for kriging model. Statistics of Models performances are available in Error! Reference source not 18 found.. On training set, MAD and RMSE values are clearly better for kriging (2.1°C and 2.9°C respectively) than for MLR (3.0°C and 4.0°C respectively). Moreover, the percentage of predicted values at a precision level

1 2

 1 Fig. 4 that all these samples are located relatively far from the rest of training data. That points the most valuable

6

  On prediction set, models performance are also quite better for kriging model than for MLR. Actually,7kriging model provides a RMSE value of 2.0°C and predicts 94% of test samples at a precision level under the 8 confidence interval of the standard measure, whereas the corresponding values are respectively 2.5°C and 75% 9 for MLR model.

Fig. 9 :4

 9 Fig. 9: Parity plots of kriging model for prediction of diesel cloud point; a) on training set; b) on

Table 1 : Example of covariance models 8

 1

	Type of covariance	Expression
	Exponential

Table 3 : Properties of the seven different VGO used to produce diesel samples 1

 3

	VGO	Origin	Density at	Nitrogen	Sulphur	IBF -FBF (°C)	Kinematic	Kinematic viscosity
			15°C				viscosity	at 100°C
			(g/cm 3 )	(wt	(wt%)		at 70°C
				ppm)
	B	Gansu	0,8974	1190,00 1,0409	348,00 -584,00	82,25	12,55
		(Chine)
	C	Iranian Saniya	0,9375	1300,00 2,8743	321,90 -615,10	124,74	11,11
	G	Husky	0,9580	1080,00 3,0300	215,00 -565,00	126,70	25,20
	K	SR	0,9346	1755,00 2,2375	337,70 -632,50	25,24	9,73
	M	Hoil+Oural	0,9314	2160,00 1,0777	292,00 -606,20	81,87	8,06
	N	Etats-Unis	0,9208	1160,00 0,2974	158,20 -581,40	32,07	4,69
	P	Mix of	0,9284	1395,00 1,8921	316,00 -624,00	133,24	11,93
		Arabian
		Light/Basrah
		(Irak) 85/15
	2
	3
	4

Table 4 : Measured and evaluated properties of feedstock and diesel samples 5

 4

	Property	Petroleum cut	Standard Methods	References
	Cloud point (CP)	diesel	NF ISO 23015	(NF EN ISO 23015
				1994)
	Simulated distillation	diesel	ASTM D2887	(ASTM D2887 2016)
	(ܶ , ݅ =
	5, 10, 15, … , 90, 95ሻ
	6
	For the prediction of diesel CP, samples were divided into training (40 samples) and prediction (16

Table 6 : Intrinsic parameters estimates in learning step for interpolation methods 10

 6

	Interpolation method	Intrinsic	Data set 1	Data set 2	Data set 3
		parameters
	Splines	Order (݉)	2	2	3
		Smoothing (ߣ)	0.17	4e-4	3e-6
	Kriging	Nugget	1.12	0.70	2e-3
		ߠ ଵ	11.50	1.24	0.45
		ߠ ଶ	10.43	1.21	0.50
	11
	3.1

.1 Data set 1 (affine function) 12

  Statistics for models performance in simulated situations are given in Table7. For data set 1, MLR

Table 7 : Statistics for Models performances in simulated situations 7

 7

	Data	Model	MRD	RMSEP
		Regression	0.13	0.24
	Data set 1	Kriging	0.07	0.40
		Spline	0.10	0.23
		Regression	7.80	2.52
	Data set 2	Kriging	3.83	1.88
		Spline	3.60	1.69
		Regression	0.40	0.32
	Data set 3	Kriging	0.10	0.18
		Spline	0.22	0.29

Table 8 : Statistics for Models performances in prediction of diesel cloud point 18

 8

	Type of	Model	MAD	RMSEP	% points	% points
	Validation		(°C)	(°C)	predicted at	predicted at
					± CI	± 2CI
	Leave-one-out	MLR	3.0	4.0	60	82
		Kriging	2.1	2.9	70	90
	Prediction	MLR	2.0	2.5	75	94
		Kriging	1.8	2.0	94	94

where ߳ denotes the measurement error related to the ݅ ௧ observation.

The simplest approach of spline so called "Interpolating spline" consists on assigning a local estimate ,ݑሺݖ ݒሻ

model which is generally based on a stochastic approach. Let ܣ denotes the study area. Thus, it is assumed that

[START_REF] Wahba | Some New mathematical -Methods for Variational Objective Analysis Using Splines and Cross Validation[END_REF]

Same statistics were evaluated on prediction set

ACS Paragon Plus Environment

Energy & Fuels

and MacOS. Various performing packages are available and allow to implement models from personal data. The 9 package "stats" contains function that allow fitting linear regression models and statistical calculations. The 10 package "DiceKriging" was used to implement kriging models (Olivier Roustant 2015). The package "fields" 11 contains a function called "Tps" that enables to develop smoothing spline models.

Appendix A : Estimation of regression parameters

In statistics, an estimator is a rule to calculate an estimate of a given quantity based on observed data [START_REF] Azaïs | Le modèle linéaire par l'exemple: Régression, analyse de la variance et plans d'expérience illustrés avec R et SAS[END_REF]. Let ߠ denotes an estimator of ߠ. It is generally obtained by solving an optimization problem such as:

where ܬ is an objective function that depends on the available observed data. There are various types of estimator according to the analytic expression of the objective function. The most currently used is the least squares (LS)

estimator which is defined as follows:

The quality of an estimator is evaluated by two statistical characteristics:

• Its bias

where ܧሾܺሿ and ܸሾܺሿ denotes respectively the expected value and the variance of the random variable ܺ. An estimator is considered as an optimal one whether it is an unbiased with a minimal variance [START_REF] Azaïs | Le modèle linéaire par l'exemple: Régression, analyse de la variance et plans d'expérience illustrés avec R et SAS[END_REF]. Under the assumptions that modelling errors are normally independently identically distributed (i. i. d.), it has been shown that the LS-estimator presents optimal characteristics (Antoniadis, Berruyer and Carmona 1992;

Azaïs and Bardet 2012).

Regression models may be divided into two groups: the generalized linear models if ݂ can be expressed as a linear combination of its parameters (affine, polynomial, quadratic, etc.); the nonlinear models (rational fraction, exponential, logarithmic, etc.)

Appendix B : multivariate normal distribution

A gaussian vector is a numeric vector such that any linear combination of his components is a normally distributed random variable. A multivariate normal distribution is entirely characterized by its mean vector ݉ and its variance-covariance matrix Σ. Let ሺܺ ଵ , … , ܺ ௗ ሻ denotes a Gaussian vector. The probability density is defined as follows:

Appendix C : General trend for universal kriging

The ݂ that appear in the definition of universal kriging predictor (paragraph Error! Reference source not found.) are polynomials of degree lower or equal to ݉. If ݔ is a ݀dimensional vector, they are denied as follows [START_REF] Journel | When do you need a trend model in Kriging ?[END_REF]: In statistics, a likelihood function is a function of a statistical model parameters given data. Likelihood functions play a key role in statistics inference. The basic idea of likelihood function consists on considering models parameters as a set of random variable. It is then defined as the probability density function of parameters given the observed data. The maximum likelihood estimate is the set of values that maximizes the likelihood function [START_REF] Myung | Tutorial on maximum likelihood estimation[END_REF].

Appendix E : Nugget effect

In the early development of geostatistics, the term 'nugget effect' was coined for the apparent discontinuity at the beginning of many semi-variogram graphs. This name was chosen to reflect the large differences found between neighboring samples in 'nuggety' mineralization such as Wits gold reefs [START_REF] Goovaerts | Geostatics for natural Resources evaluation[END_REF]. It may be defined as follows:

where ܥ is the recognized covariance model by plotting semi-variogram.


Baptiste Magnier
email: baptiste.magnier@mines-ales.fr

Philippe Montesinos
email: philippe.montesinos@mines-ales.fr

Oriented Half Gaussian Kernels and Anisotropic Diffusion

Keywords:

Nonlinear PDEs (partial differential equations) offer a convenient formal framework for image regularization and are at the origin of several efficient algorithms. In this paper, we present a new approach which is based (i) on a set of half Gaussian kernel filters, and (ii) a nonlinear anisotropic PDE diffusion. On one hand, half Gaussian kernels provide oriented filters whose flexibility enables to detect edges with great accuracy. On the other hand, a nonlinear anisotropic diffusion scheme offers a means to smooth images while preserving fine structures or details, e.g. lines, corners and junctions. Based on the calculus of the gradient magnitude and two diffusion directions, we construct a diffusion control function able to achieve precise image regularization. Some quantified experimental results compared to existing PDEs approaches and a discussion about the parameterizing of the method are presented.

A Framework of Anisotropic Diffusion with PDE

Obtain regularized versions of noisy, corrupted or degraded images caused for example by compression artifacts is a difficult task in image processing. However, preserving significant internal structures is a field that has largely benefited from techniques of Partial Differential Equations (PDE) (Aubert and [START_REF] Kornprobst | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF]Magnier and Montesinos, 2013). PDEs belong to one of the most important part of mathematical analysis and are closely related to the physical world. In this context, images are considered as evolving functions of time and a regularized image can be seen as a version of the original image at a special stage. Thereby, PDEs methods smooth locally the image following one or several directions which are different in each point of the image. In this paper, let us note I : Ω → R, (Ω ⊂ R 2 ) a grey level image with I(x, y) corresponding to the pixel intensity of coordinates (x, y). Considering I 0 the original image, the general evolution model can be formally written in the following form: ∂I ∂t (x, y,t) = F (I(x, y,t)) I(x, y, 0) = I 0 (x, y)

(1)

where F is a given image processing algorithm, preserving edges having high gradient. F represents a function of the original image I 0 and its first and second order spatial derivatives [START_REF] Caselles | Introduction to the special issue on partial differential equations and geometrydriven diffusion in image processing and analysis[END_REF].

It should be noted that Koenderink [START_REF] Koenderink | The structure of images[END_REF] was the first to underline the equivalence between the convolution with a Gaussian kernel of standard deviation √ 2t and the solution of the PDE describing the heat diffusion, at a time t. This smoothing process, called isotropic diffusion, is known to smooth noise but blur edges, leading to loose image structures. In order to regularize images by controlling the diffusion, Perona and Malik [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF]) have proposed a model described by the following equation:

∂I ∂t (x, y,t) = div (g ( ∇I ) • ∇I ) (2)
where div represents the divergence operator and g(s) : [0, +∞[→]0, +∞[ a decreasing function satisfying g(0) = 1 and g(+∞) = 0, this function could be chosen as

g ( ∇I ) = e - ∇I K

2

, with K ∈ R a constant that can be assimilated to a gradient threshold or a diffusion barrier. The decomposition of the eq. 2 with the second derivatives of I in orthogonal directions (ξ ⊥ η) respectively in the edge direction called ξ and in the gradient direction labelled η = ∇I ∇I enables to understand the diffusion behavior [START_REF] Kornprobst | Nonlinear operators in image restoration[END_REF]:

∂I ∂t (x, y,t) = c ξ • I ξξ + c η • I ηη (3)
where (I ξξ , I ηη ) = ∂ 2 I ∂ξ 2 , ∂ 2 I ∂η 2 , c ξ and c η are coefficients tuning the diffusion (diagrammed in Fig. 1(a)).

When c ξξ = c ηη = 1, the eq. 3 is equivalent to the heat equation [START_REF] Koenderink | The structure of images[END_REF]. Choosing c ξ = g ( ∇I ), a gradient function and c η = g ( ∇I ) + ∇I • g ( ∇I ), the diffusion process described in eq. 3 can be interpreted as two directional heat flows 1 with different diffusion intensities in the η and ξ directions to preserve discontinuities:

• Inside homogeneous regions, the gradient magnitude ∇I is small and the diffusion becomes isotropic.

• On edges, the diffusion becomes anisotropic, being attenuated by the function g, and is inhibited when the two coefficients (c ξξ , c ηη ) tend to zero.

Diffusion control is done with finite differences so that many contours of details are preserved. However, within images corrupted by a heavy noise, generally, this noise is not totally removed because the diffusion process is inhibited.

Gaussian filtering for gradient estimation has been used in a number of works to elaborate the model presented in eq. 3 less sensitive to noise and more stable. We can mention here the approach of Alvarez et al. [START_REF] Alvarez | Image selective smoothing and edge detection by nonlinear diffusion, ii[END_REF] which induces for each pixel either an adaptive unidirectional tangential diffusion I ξξ at level of edges or an efficient isotropic smoothing for noise removal inside homogeneous regions. Nevertheless, this smoothing model does not allow a progressive diffusion in the gradient direction η because it depends on two diffusion barriers. Consequently, in the presence of a high noise, even in homogeneous regions, this diffusion scheme behaves like the Mean Curvature Motion [START_REF] Catté | Image selective smoothing and edge detection by nonlinear diffusion[END_REF] (MCM) method which consists in performing the diffusion only along the tangential direction ξ or along isophote lines. Although the MCM scheme regularizes the image in edge directions, this approach tends to round corners after a certain number of iterations and can create stripes inside noisy homogeneous regions.

Instead of considering only the gradient magnitude to drive the diffusion, tensorial approaches [START_REF] Weickert | Coherence-enhancing diffusion filtering[END_REF][START_REF] Tschumperlé | Vector-valued image regularization with pdes: A common framework for different applications[END_REF][START_REF] Tschumperlé | Fast anisotropic smoothing of multi-valued images using curvature-preserving PDE's[END_REF] contribute to another image diffusion formalism. From a structure tensor J ρ = G ρ * ∇I σ ∇I T σ , where G ρ denotes a Gaussian kernel of standard deviation ρ, authors of [START_REF] Weickert | Coherence-enhancing diffusion filtering[END_REF][START_REF] Tschumperlé | Vector-valued image regularization with pdes: A common framework for different applications[END_REF][START_REF] Tschumperlé | Fast anisotropic smoothing of multi-valued images using curvature-preserving PDE's[END_REF] elaborate a tensor field which specifies the local smoothing geometry defined from the spectral elements of the structure tensor. Then, us-

1 Note that if ∇I > K √ 2
, then c η <0 and the diffusion equation behaves locally like an inverse diffusion equation which is an unstable process enhancing features. ing the divergence [START_REF] Weickert | Coherence-enhancing diffusion filtering[END_REF] or the trace [START_REF] Tschumperlé | Vector-valued image regularization with pdes: A common framework for different applications[END_REF], the smoothing along a contour in inversely proportional to the contour strength in the direction of the eigenvector associated to the higher eigenvalue. Inside homogeneous regions, eigenvalues are close to zero and the diffusion becomes isotropic.

As demonstrated in [START_REF] Tschumperlé | Vector-valued image regularization with pdes: A common framework for different applications[END_REF], trace based PDE is best suited to understand the local smoothing geometry behavior and these diffusion scheme ensure coherence smoothing directions but the Gaussian behavior on curved structures or corners results in a "mean curvature flow effect" leading to round small structures or corners. In order to compensate this drawback, the author of [START_REF] Tschumperlé | Fast anisotropic smoothing of multi-valued images using curvature-preserving PDE's[END_REF] proposed a curvature-preserving smoothing PDE that diffuses the image I along a field of vectors w issued by the eigenvectors of J ρ . Despite the fact that the author of this method [START_REF] Tschumperlé | Fast anisotropic smoothing of multi-valued images using curvature-preserving PDE's[END_REF] has demonstrated that it better preserves corners and small details in the image, as the other tensorial approaches, when the anisotropic coefficient is too large, the diffusion of a high noise brings a fiber effect in homogeneous regions. To avoid this undesired diffusion effect, it is preferable to use a higher standard deviation of the Gaussian σ, however this leads to delocate even so the corners, diffuse small objects and also blur edges.

In this paper, we propose a new PDE scheme that regularizes images considering two contour directions. This diffusion process correctly preserves corners as well as small objects and becomes isotropic inside homogeneous regions without generating undesired fiber effect of artifacts. Thanks to a rotating Gaussian derivative half-filter, we extract a gradient amplitude and determine two edge directions.

As illustrated in Fig. 1(b), if the gradient magnitude is small, the two major smoothing directions are straightened for an opposite alignment (i.e. 180 • ), becoming thereby either tangential to the edge or creating an isotropic diffusion inside flat regions. Then we apply an anisotropic diffusion on each pixel using new control functions adapted to our objectives of image regularization.

Oriented Half Gaussian Kernels

Steerable filters [START_REF] Freeman | The design and use of steerable filters[END_REF][START_REF] Jacob | Design of steerable filters for feature detection using canny-like criteria[END_REF] or anisotropic edge detectors [START_REF] Perona | Steerable-scalable kernels for edge detection and junction analysis[END_REF] perform well in detecting large linear structures (see kernels in Fig. 2 (a) and (b)). Close to corners however, the gradient magnitude decreases as the edge information under the scope of the filter decreases. Consequently, the robustness to noise concerning small objects becomes very weak.

A simple solution to bypass these effects would be to consider paths crossing each pixel in several directions. Wedge steerable filters introduced by Simoncelli and Farid [START_REF] Simoncelli | Steerable wedge filters for local orientation analysis[END_REF] are composed of asymmetric masks providing orientation of edges in different directions issued from a pixel. Akin to oriented histograms, the saliency of the gradient measure is calculated at each of discretized orientations (a wedge filter is presented in in Fig. 2 (d)). An advantage of these oriented filters is that they allow a characterization of junctions [START_REF] Mühlich | Design and implementation of multisteerable matched filters[END_REF]. In [START_REF] Michelet | Estimating local multiple orientations[END_REF] is presented also an asymmetric operator based on a sliding rectangular window where the orientation is defined as being the angle that corresponds to the maximum homogeneity i.e. the minimum variance. Unlike the Gaussian function, which is an optimal solution for the Canny criteria [START_REF] Canny | A computational approach to edge detection[END_REF], in the direction of the edges, these oriented filters have a constant amplitude on almost the whole extent of the mask.

The idea developed in [START_REF] Montesinos | A New Perceptual Edge Detector in Color Images[END_REF] was to "cut" the derivative (and smoothing) Gaussian kernel in two parts: a first part along an ini-

Q (x, y, θ) = I θ * C 1 • H (-y) • x • e - x 2 2λ 2 + y 2 2µ 2 (4)
where I θ corresponds to a rotated image2 of orientation θ, C 1 is a normalization coefficient, (x, y) are pixel coordinates, and (µ, λ) the standard deviations of the anisotropic Gaussian filter. Since we only require the causal part of this filter along Y axis, we simply "cut" the smoothing kernel by the middle, in an operation that corresponds to the Heaviside function H. Q (x, y, θ) represents the slope of a line derived from a pixel in the perpendicular direction to θ.

To obtain gradient magnitude measure ∇I and its associated direction η on each pixel P, we first compute the global extrema of the function Q (x, y, θ), with θ 1 and θ 2 . θ 1 and θ 2 define a curve crossing the pixel (an incoming and outgoing direction). Two of these global extrema are combined to obtain ∇I :

         ∇I = max θ∈[0,360[ Q (x, y, θ) -min θ∈[0,360[ Q (x, y, θ) θ 1 = arg max θ∈[0,360[ (Q (x, y, θ)) θ 2 = arg min θ∈[0,360[ (Q (x, y, θ))
(5) Once ∇I , θ 1 and θ 2 have been obtained, the edges can be easily extracted by computing local maxima of ∇I in the direction of the angle η = (θ 1 + θ 2 )/2 followed by an hysteresis threshold (see [START_REF] Montesinos | A New Perceptual Edge Detector in Color Images[END_REF] for further details). In this paper, we are solely interested in the gradient magnitude and the two directions (θ 1 , θ 2 ) which are improved (see next section) to be used in our diffusion scheme discussed in Section 3.2.

Due to their adjustable lengths, rotating filters enable computing two precise diffusion orientations in the edge directions, even at high noise levels. More details about the effect of noise can be found in (Magnier et al., 2011a) where the authors have evaluated this edge detector as a function of noise level, and, a comparison with other approaches [START_REF] Perona | Steerable-scalable kernels for edge detection and junction analysis[END_REF] shows the efficiency of this method. Note that these kernels have been used in several diffusion schemes in image regularization [START_REF] Magnier | Texture Removal in Color Images by Anisotropic Diffusion[END_REF][START_REF] Magnier | A new region-based pde for perceptual image restoration[END_REF]Magnier et al., 2013a;[START_REF] Magnier | Perceptual color image smoothing via a new region-based pde scheme[END_REF]Magnier and Montesinos, 2013).

3 Smoothing in Two Edge Directions

Two Improved Edge Directions in a Diffusion Sceme

As detailed in Section 1, PDE-based image regularization techniques using gradient intensities or tensorial diffusion smooth the image either in the orthogonal directions (ξ, η), or in the directions provided by the eigenvectors of the structure tensor. However, both approaches do not take into account the two actual directions of edges, for example at a corner level.

For removing texture and preserving edges, issued by an edge classifier, the original idea developed in [START_REF] Magnier | Texture Removal in Color Images by Anisotropic Diffusion[END_REF] was to smooth the image in the two contour directions called (ξ 1 , ξ 2 ):

∂I ∂t (x, y,t) = I ξ 1 ξ 2 = ∂ 2 I ∂ξ 1 ∂ξ 2
. As textures are annihilated, this approach is not adapted for image regularization, especially because the diffusion is not controlled.

In [START_REF] Magnier | Perceptual color image smoothing via a new region-based pde scheme[END_REF], a pixel classification determines roughly if a pixel belongs to a homogenous region or an edge, then, authors have developed a new diffusion method. Inside edge regions, a function of the gradient magnitude issued by half Gaussian kernels (eq. 5) and also of the angle between the two diffusion directions (θ 1 , θ 2 ) called β = abs(θ 1θ 2 ) (Fig. 4(a)) drives the diffusion process:

∂I ∂t (x, y,t) = e -∇I K 1 2 + e -(π-β) (π•K 2 ) 2 2 • ∂ 2 I ∂θ 1 ∂θ 2 , (6)
with K i,i∈{1,2} ∈ ]0; 1]. This diffusion process combines isotropic and anisotropic diffusion, while preserving the edges and corners of different objects in highly noisy images. Nevertheless, instead of enhance textures, authors use of the heat equation, smoothing them isotropically. Moreover, as the anisotropic diffusion process is applied only at position of edges, it creates undesired lineaments near edges caused by the directional diffusion (θ 1 , θ 2 ).

The anisotropic edge detector based on half kernels informs (eq. 5) both on the importance of the contours (gradient value) and its directions (θ 1 , θ 2 ). Our regularization method uses these two orientations in the smoothing process in function of the angle β between θ 1 and θ 2 . However, inside homogeneous regions, due to the level lines on the image surface, θ 1 and θ 2 are not opposite directions (i.e. β = π, see Fig. 6(c)), thus the diffusion does not perform isotropically. Here, we propose to align the orientations diffusion when the gradient value is low (i.e. θ 1 = θ 2 + π). To achieve this, a new angle called γ between the two edges directions is defined as:

γ = min(β + π • e - ∇I a 2 -π • e -1 a 2 , π) (7)
The more ∇I is close to zero, the more γ is close to π (plotted in Fig. 4(a) bottom). The parameter a ∈ ]0, 1] steers the straightening, as shown in Fig. 4, the more a is growing, the more γ is close to π for each pixel. Now, as illustrated in Fig. 4(a), we consider the two new diffusion directions ρ 1 and ρ 2 in our method (see Section 3.2), estimated modulo 2π in function of the following table :

θ 1 > θ 2 θ 1 < θ 2 ρ 1 = θ 1 + γ-β 2 ρ 1 = θ 1 -γ-β 2 ρ 2 = θ 2 -γ-β 2 ρ 2 = θ 2 + γ-β 2

An Approach Preserving Contours and Smoothing Regions

Our algorithm enables a smoothing in the improved directions of the contours, thus preserving edges and details (I ρ 1 ρ 2 term) while diffusing also in the direction of η for edges having a low gradient or inside homogeneous regions (I ηη term). Furthermore, these three directions smoothing terms have to be controlled in order to preserve image contours and not to create undesired artifacts or fiber effect elsewhere. In The more the image must be regularized, the more the a parameter must be large (e.g. a = 0.2 for a heavy noise, see bottom). The more the pixel value is dark, the more the ange is acute.

this respect, we adapt a new PDE developed in (Magnier and Montesinos, 2013), involving the gradient magnitude (eq. 5) and the γ angle (eq. 7) driving both the diffusion terms I ρ 1 ρ 2 and I ηη :

               ∂I ∂t (x, y,t) = f k • I ρ 1 ρ 2 + f h • I ηη f k = e - ∇I k 2 + e - (π-γ) (π•k) 2 2 f h = e - ∇I h

2

, with (k, h) ∈ ]0; 1] 2 (8) The f k function ensures the diffusion preserving edges and corners whereas the f h function enables a permanent smoothing in the gradient direction for noisy homogeneous regions. One one hand, the more the h value is close to 1, the more edges are blurred, one the other hand the more the k value is close to 1, the more the diffusion process is important. Contrary to [START_REF] Alvarez | Image selective smoothing and edge detection by nonlinear diffusion, ii[END_REF], these control functions are not threshold functions but continuous functions (Fig. 5). Thus, the diffusion is never only in the (ρ 1 , ρ 2 ) directions. In case of a small gradient and a γ angle close to π, the considered pixel will be widely diffused (see Fig. 5(a)). If the gradient absolute value is large and the γ angle is small, smoothing is weak and operates mainly along these two orientations (ρ 1 , ρ 2 ).

Experimental Results and Analysis

In this section, we present several results of our regularization method compared to different other PDEs approaches. For each result, presented below, we detail the parameters used either for our algorithm, or for other methods. Also, we show a SSIM [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF]) evaluation of noisy images as a function of the number of iterations using different parameters, which permits us to discuss about the choice of the best parameters couple (k, h). Note that, in order to obtain precise diffusion directions (θ 1 , θ 2 ), we use a discretization angle of ∆θ = π 90 = 2 • and the standard deviation of the half anisotropic Gaussian are µ = 5, λ = 1 for the gradient extraction (eq. 5).

Firstly, the image in Fig. 6(a) is corrupted by a low acquisition noise because it is a scan of the original Lena from Playboy. After 10 iterations, our method removes this noise, and plain regions become totally homogeneous whereas small objects are well preserved (see details in Fig. 6(b)). Diffusion parameters are (k, h)=(0.3, 0.1) and a = 0.1 because the noise is low in this image.

The second image shown in Fig. 8(b) is a natural image contaminated by a Gaussian noise (σ = 10). We aim to regularize this picture and preserve edges as far as possible. As presented in Fig. 8(g),(h) and (i), comparing the absolute error between the original image and the regularized image, our algorithm preserves better edges than tensorial results [START_REF] Tschumperlé | Vector-valued image regularization with pdes: A common framework for different applications[END_REF][START_REF] Tschumperlé | Fast anisotropic smoothing of multi-valued images using curvature-preserving PDE's[END_REF]. To obtain a better visualization, note that the absolute error images are corrected following a curve on the image histogram, as presented in Fig. 8(c). This visualization process is the same for each absolute error image presented in this paper.

The third noisy image presented in Fig. contains a random Gaussian noise of standard deviation σ = 20. Due to this high noise and the texture, this image is particularly difficult to regularize correctly while preserving the thin texture. For each filter, we choose the parameters that gives the best results. In order to obtain comparative results, we choose the same width (i.e. standard deviation) of the Gaussian kernel for approaches using this function (i. e. σ = µ = 1 for [START_REF] Alvarez | Image selective smoothing and edge detection by nonlinear diffusion, ii[END_REF][START_REF] Weickert | Coherence-enhancing diffusion filtering[END_REF][START_REF] Tschumperlé | Vector-valued image regularization with pdes: A common framework for different applications[END_REF][START_REF] Tschumperlé | Fast anisotropic smoothing of multi-valued images using curvature-preserving PDE's[END_REF]). We compare our result with the MCM [START_REF] Catté | Image selective smoothing and edge detection by nonlinear diffusion[END_REF], the methods of Perona-Malik [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF], Alvarez et al. [START_REF] Alvarez | Image selective smoothing and edge detection by nonlinear diffusion, ii[END_REF], tensorial driven diffusion [START_REF] Weickert | Coherence-enhancing diffusion filtering[END_REF][START_REF] Tschumperlé | Vector-valued image regularization with pdes: A common framework for different applications[END_REF][START_REF] Tschumperlé | Fast anisotropic smoothing of multi-valued images using curvature-preserving PDE's[END_REF] and Magnier et al. [START_REF] Magnier | Perceptual color image smoothing via a new region-based pde scheme[END_REF].

It is easy to remark that MCM and Perona-Malik models do not remove correctly the noise. Algorithm of Alvarez et al. loses the texture and creates artifacts at position of edges. Tensors [START_REF] Weickert | Coherence-enhancing diffusion filtering[END_REF][START_REF] Tschumperlé | Fast anisotropic smoothing of multi-valued images using curvature-preserving PDE's[END_REF] diffusion creates a fiber effect in homogeneous regions due to the high noise. Tensorial result of Tschumperlé [START_REF] Tschumperlé | Vector-valued image regularization with pdes: A common framework for different applications[END_REF] gives a good result even if this method is known to distort corners but creates a graininess effect inside homogeneous regions (see details in Fig. 9(l)). Our algorithm restores correctly edges (Fig. 9(i)), does not create undesirable fiber effect inside homogeneous regions and enables a regularization of stripes textures with the use of thin half kernels (details in Fig. 9(m)).

We have tested different values of the couple (k, h) and when k < h, the diffusion results in blurring edges. When k 0.2, flat regions become totally homogeneous but the smoothing process is so weak on edges that they are not retored. Fig. 7 enables a better visualization for the choices of the couple (k, h) in function of the SSIM evolution. It determines that for (k, h) = (0.3, 0.1), the results are the best with this set of parameters, only the a parameter (Fig. 4) and the standard deviations of the half kernels must be set in function of the image type. Generally, the choice of our half kernel filters parameters (µ, λ) has to be done in function of the noise level. In order to preserve small objects, we can choose for the length of our filter µ = 5. However, the width which corresponds to the derivation filter depends on the noise level and the image size. typically, we can choose λ = 1, it enables to keep precisely edges but if the noise is higher, we can choose a larger value.

Conclusion

In this paper, we have presented a new regularization technique based on half Gaussian filtering followed by an anisotropic PDE diffusion. This method complements previous works by calculating new directions of diffusion for a more accurate regularization of images. We have shown that the fine assessment of the diffusion direction is a key factor, and we have proposed a set of control parameters which provide very good results in restoration of noisy images. Eventually, as these diffusion directions preserve precisely corners and small objects in images, future works will focus on the extraction of corners only using the improved edges directions formula detailed in this paper.

Figure 1 :

 1 Figure1: Diagrams of edge diffusion. (a) An image contour and its moving vector basis (ξ, η) and diffusion representation with ellipsoids. The more the gradient is high, the more the ellipse is elongated. Gradient and tangential direction denoted (η, ξ) and diffusion representation with ellipsoids, note that ellipsoids are not always oriented in the ξ direction using tensorial methods. (b) Our desired diffusion representation with half ellipsoids, the more the edge is sharped and the angle is acute, the more the half ellipses are thin.

  Figure 2: Different 2D derivative Gaussian kernels and a wedge filter

Figure 5 :

 5 Figure 5: Control functions of the diffusion process.

Figure 4 :

 4 Figure 4: Different angles representation, the straightening up function and γ angle representation for corrupted images in Figs. 8(b) and 9(b) in function of different values of a.The more the image must be regularized, the more the a parameter must be large (e.g. a = 0.2 for a heavy noise, see bottom). The more the pixel value is dark, the more the ange is acute.

Figure 6 :

 6 Figure 6: Our algorithm applied on a real image polluted by an acquisition noise.

Figure 8 :

 8 Figure 8: Image regularization and absolute error. Absolute error are negative images.

As explained in[START_REF] Montesinos | A New Perceptual Edge Detector in Color Images[END_REF], the image is oriented instead of the filter (like the oriented filter presented in[START_REF] Michelet | Estimating local multiple orientations[END_REF]) so as to decrease algorithmic complexity and to allow use of a recursive Gaussian filter[START_REF] Deriche | Recursively implementing the gaussian and its derivatives[END_REF]. As a matter of fact, for implementation purpose we replace the filtering of the initial image I 0 by a filter oriented along the varying direction θ with the filtering of an image I θ , rotated by an angle -θ, by the constant filter of orientation 0. This last operation is described by eq. (1) and is completely equivalent to a rotated filtering.