La production de cette thèse est avant tout le fruit d'un travail d'équipe basé sur la confiance entre personnes. Pour cela je remercie Anne-Marie HABRAKEN (

Je remercie mes deux supers garçons Clarence-Félix et Simon-Maxence, dont j'ai la chance de les voir grandir au point de devenir de vrais dictionnaires Franco-allemands vivants. Aujourd'hui vous ne comprenez pas tout à fait ce que je faisais en Allemagne, en Belgique et en France durant ces trois années mais dans quelques années vous aurez la réponse en lisant ce manuscrit. Alors encore un peu de patience mes garçons. Ce travail est aussi pour que vous sachiez qu'après être devenu athlète de haut niveau, éducateur dans un centre de jeunesse, musicien, entrepreneur, ingénieur, j'étais apprenti scientifique.

En espérant n'avoir oublié personne je tiens à remercier une personne sans qui cette thèse n'aurait jamais vu le jour, c'est ma femme, la courageuse, grande de taille et de coeur, Silvia FANSI.

Etat d'esprit adopté durant la période de thèse: Le sentiment d'isolement que j'ai parfois ressenti lorsque j'étais loin de ma famille a été un atout efficace pour me concentrer sur l'essentiel, à savoir, ma femme, mes enfants, ma famille en France et mon travail de Thèse. Ce travail dont j'ai investi beaucoup de temps et d'énergie est en fait le point de départ et la base d'une réflexion pour que dans le future mes enfants et ma famille puissent croire en des choses qu'ils veulent réaliser et de s'engager jusqu'au bout, en travaillant très dur pour cela. Je suis persuadé que l'intelligence ou le savoir n'a ni odeur, ni couleur et encore moins de limite. Il n'existe pas une intelligence mais plusieurs. Il existe une expression qui résume cet état d'esprit et que les africains aiment fréquemment utiliser dans des situations particulièrement difficiles "On va faire comment?".

~I.3 ~

I.I.3 Dual-Phase in automotive industry

As already mentioned, the selected DP steels in the present study are the widely advanced high-strength steels used in the automotive industry confirms that, DP steels usually consist of a ferrite matrix containing a hard martensite second phase in the form of islands. Increasing the volume fraction of hard second phases generally increases the strength. DP steels are obtained by controlled cooling from the austenite phase (in hot-rolled products) or from the two-phase ferrite plus austenite phase, to transform some austenite to ferrite before a rapid cooling transforms the remaining austenite to martensite. According to [PAR 2012], roll forming, as a metal fabricating process, is used to add both strength and rigidity to manufactured DP steels. Figure I.4 illustrates the application of various DP steels in automotive industry such as bumper beam, A-frame reinforcement, roof bow, B-pillar reinforcement, rear side member, front floor cross member, floor side reinforcement, front side members and Body-in-White (BIW) structural parts.

~I.5 ~

I.2 Objectives of the thesis

DP steels modelling has been since 1980 a huge interest for car and steel makers. The control of their mechanical properties depends on their formability characteristic which influences the manufacturing quality. Regarding the constant increasing complexity of the technology employed to reach the industrial requirements, the needs of experiment and analytical models are lead by the development of numerical simulations. This technology, born forty years ago is in constant improvement years after years thanks to powerful computers and broad research on the Finite Element method and on constitutive laws based on phenomenological or multiscale approaches. Nowadays, the accuracy of the simulations results is gaining credibility compared to the experimental test. The numerical simulation reduces dozen of experiments, time framework and manufacturing costs. On top of that, an accurate simulation needs realistic boundary conditions, good sketch of the sample and tooling geometries and the most important a reliable material behaviour model. Regarding the last requirement, the needs for ArcelorMittal is to predict with good accuracy the rupture in forming processes taking into account the edge effects, the residual stress, the ~I.6 ~ The main goal of this work is to correlate the experimental results on notched samples and the model predictions for better understanding of the DP steel ductile fracture. To attempt this aim, the GTNB model has been adapted as "User-defined Material model subroutine" (VUMAT) in the Abaqus/explicit FE code [ABA 2011]. The model has been enriched by adding a coalescence model, a recent void nucleation and growth laws integrating the back stress variable [LAN 2011] and a fracture initiation criterion. These enhancements have been done based on high resolution X-ray tomography observations and measurements. In order to accurately correlate the finite element predictions with the experiments, a precise postprocessing method has been developed taking into account identical results extraction between experiments and simulations. The numerical void volume fraction definition is the same as the one used in the test related to the number of cavities and their size in a fixed volume. Along with the attempt of integrating the back stress variable in the new nucleation law, a discussion has been opened regarding the triaxiality definition. The second objective is to be able to rank the DP steel grades chosen in this study. Finally, the experimental validation is being further extended to other sample geometries in flat sheet steels, as well as industrial application such as cross-die drawing test.

~I.7 ~

I.3 Contents

The aim of the second chapter is to give a general presentation of the Dual-Phase (DP) steels. It will show typical DP steels microstructures, mechanical properties and formability characteristics. The third chapter presents how the mechanical behaviour of DP steels can be predicted by plasticity, damage or fracture modelling. The experiments carried out on in-situ tensile notched axisymmetric specimens of DP steels [LAN 2011] revealed a strong dependency between the density of voids, the back stress, and the triaxiality. The fourth chapter summarizes the contributions brought to the damage and fracture modeling during this thesis. These contributions are based on the recent experimental results of [LAN 2011] on in-situ tensile notched axisymmetric specimens of DP steels, and the onedimensional metallurgical models inferred therein. These contributions concern the void nucleation and void growth law, the coalescence model and an additional fracture initiation criterion, which were extended to the 3D incremental case and integrated in the proposed model (further called GTNBF). The parameter identification of the whole set of parameters of the GTNBF model for three DP steels is described in chapter five. Chapter six deals with the finite element implementation of proposed model in the commercial code Abaqus/Explicit. The implementation follows the path previously adopted by Ben Bettaieb [BEN 2012] in the finite element code Lagamine. The proposed implementation remains compatible with the two finite element codes, with the Abaqus routine "VUMAT" being used as an interface. The developed algorithm and code is numerically validated. The potentialities and limitations of the model are contained in Chapter seven where the material parameters sensitivity is carefully checked. In front of that, the post processing and the element size influence on the current model are carried out. Chapter eight illustrates the industrial applications of the GTNBF model. Finally, the conclusions are summarized in chapter nine along with some potential future work.

~II.2 ~

The aim of the current chapter is to give a general presentation of the Dual-Phase (DP) steels. The first part presents the typical properties of DP steels: microstructures, mechanical properties and formability characteristics. The second part of the chapter focuses on the experimental damage tests realized to evaluate the void evolution in DP steels, and gives more details of the DP steels studied. observations, the ferrite is continuous for many grades up to DP780 steel, but as volume fraction of martensite exceed 50 percent (as might be found in DP980 steel or higher strengths), the ferrite may become discontinuous.

II.1 Description of DP steels

In addition to the presence of this residual phase, the DP steels can contain negligible quantities of inclusion particles or voids generated during the production of the material (cold or hot-rolling procedures [TSI 2006] ) or after the forming process of the DP steel blank.

II.1.2 Mechanical properties

The variety of microstructures and especially the volume fraction of martensite reveal a large multitude group of tensile strength levels. Table II

II.I.3 Formability characteristics

The Dual-Phase steels present an excellent candidate for the car body structural components. These are often produced for safety-critical parts (see Figure II.4) which maintain passenger surviving space in crash events. The DP steels present a good balance of strength, formability, energy absorption and durability. Also the employability of this kind of steel provides the possibility of reducing the weight of the vehicle. DP is sometimes selected for visible body parts and closures, such as doors, hoods, front and rear rails. Other well known applications include: beams and cross members; rocker, sill, and pillar reinforcements; cowl inner and outer; crush cans; shock towers, fasteners, and wheels [TAM 2011]. Another crucial cause (not studied in this document) for the fabrication of inconsistent sheet metal part is springback i.e. the elastic strain recovery in the DP steel after the tooling is removed. Springback of sheet metal parts after forming causes deviation from the designed target shape and produces downstream quality problems and assembly difficulties as seen in ~II.9 ~

II.2.3 In-situ tensile tests

Due to the fact that the X-ray tomography is a non destructive test as mentioned previously, this technique has been coupled with an in-situ tensile test as seen in

III.1 Introduction

The previous chapter presented the DP steels as the coexistence of hard islands of Martensitic phase dispersed in a soft Ferritic matrix. This microstructure provides a good combination of strength and ductility. The mechanical behaviour until failure has been studied for more than four decades. Nowadays, the ductile failure mechanism is well known and divided in three steps before fracture: void nucleation, void growth, and void coalescence. Chapter III describes the different proposals that the scientific community has developed to model the mechanical behavior of metals such as elasto-plastic damage constitutive laws as well as rupture criteria. Of course as often as possible, examples will be dedicated to DP steels.

Many researchers have contributed to give a proper microscopic approach. Habraken [HAB 2004] presents general features of crystal plasticity models and homogenization techniques to reach macroscopic scale while [KAD 2011], [VAJ 2012], [CHO 2013] are studying specifically DP models at microscopic or macroscopic scales. Very often macro scale is the world of phenomenological approaches which can however have roots within micro scale. The mechanical behaviour of DP steels can be modeled, based on elasto-plastic theory including both scales.

The phenomenological approach of elasto-plastic behaviour will be used within this thesis. It can mainly be defined by three different assumptions: a yield function, a hardening model and a plastic flow rule.

The first hypothesis is described by an initial yield surface defined in stress space. This function also called plastic yield criterion is a mathematical description of the initial yield surface. It can be isotropic (von Mises [MIS 1928], Tresca [TRE 1868]) or anisotropic (Hill48 [HIL 1948], Barlat [BAR 2004]).The second hypothesis known as hardening model describes the evolution of the shape, the size and the position of the yield surface during the deformation. It is mainly divided in two categories: isotropic and kinematic hardening. The isotropic hardening models the expansion of the yield surface with no shape distortion while the kinematic hardening also called anisotropic hardening computes the yield surface displacement in the stress space. Shape distortion is only seldom addressed with phenomenological models.

The third hypothesis called the flow rule defines the relation between the plastic strain rate tensor and the stress tensor. A plasticity model is called associative, if the yield function is considered as a plastic potential and its derivative provides the strain rate direction. the microscopic approach, not studied within this thesis, the global macroscopic stress and strain tensors are calculated relying on a numerical simulation of each particle or grain of the DP steels [KAD 2011], [VAJ 2012], [CHO 2013]. This approach physically describes the heterogeneity of plastic strain contained in the material. However, it requires a huge quantity of data storage and CPU time. Both approaches are complementary, the microscopic approach

  Figure I.3: Strength-Elongation curve for low strength, conventional HSS and first generation AHSS steels [WAG 2006].

Figure I. 4 :

 4 Figure I.4: Automotive industry applications of the DP steels in automotive industry [PIN 2012].

  material and mechanical anisotropy. Therefore a study has been carried out in 2008 by a joint research composed of MS²F Argenco (Liège, Belgium) and ArcelorMittal-Maizières (Maizières-Lès-Metz, France) targeting the development of a formability criterion for deep drawing, rolling, bending etc... As a result of this research study conducted by Ben Bettaieb for 2 years, an extended Gurson-Tvergaad-Needleman (GTN) model [BEN 2010] has been proposed for ductile failure. The study included a literature review and state of the art of formability criteria [BEN 2011], the development of a 3D advanced GTN Formability criterion (damage modelling) [BEN 2012] implemented in the finite element code Lagamine [ZHU 1992]. The Gurson-Tvergaard-Needleman-Ben-Bettaieb model (GTNB) takes into account the void nucleation. It allows for the accurate modelling of the observations of damage initiation and growth in DP steels measured by high-resolution X-ray absorption tomography [BOU 2008]. The numerical implementation used an explicit-implicit algorithm. It is explicit for the porosity state variables and implicit for other variables (macroscopic plastic strain and the yield stress of the dense matrix, backstress, equivalent plastic strain). This choice, justified by the complexity of the porosity function and dependency of the other variables, requires small time steps to avoid any convergence problem during the computation. The success of the GTNB model application on smooth specimen motivates the joint research to stretch the model out and to implement it into a commercial finite element code used by ArcelorMittal, through a three years PhD work. At the starting point of the present research, further X-ray tomography measurements have been investigated by Landron et al. [LAN 2011] on in situ tensile notched specimens made of DP steels. The experiments revealed a strong dependency between the density of voids, the back stress, and the triaxiality for these grades. Motivated by these new experimental observations and the industrial needs, the extension of the GTNB model has been developed in cooperation with four research teams (ArcelorMittal-Maizières, INSA-Lyon, ENSAM ParisTech-Metz and University of Liège) [FAN 2013].

1 II. 1

 11 Chapter II. Dual-Phase steels Contents Chapter II. Dual-Phase steels ............................................................................ II.Description of DP steels .......................................................................... II.2 II.1.1 Microstructures............................................................................... II.2 II.1.2 Mechanical properties .................................................................... II.3 II.I.3 Formability characteristics ............................................................. II.5 II.2 Experimental Damage Investigation of DP steels .................................. II.7 II.2.1 Experimental techniques to study fracture ..................................... II.7 II.2.2 X-ray tomography principle ........................................................... II.8 II.2.3 In-situ tensile tests .......................................................................... II.9 II.3 Selected DP steels grades ...................................................................... II.10

II. 1 . 1 Microstructures

 11 Dual Phase steels belong to the larger category of Advanced High Strength Steels (AHSS) used by automakers. They provide an outstanding combination of strength and ductility as a result of their microstructure, in which hard islands of Martensitic phase (α') are dispersed in a soft Ferritic matrix (α) as shown in Figure II.1.

Figure II. 1 :

 1 Figure II.1: SEM micrographs of a typical DP steel's microstructure showing (a) ferrite (α)matrix along with banded islands of martensite (α'), (b) sub-structure within (α') phase, and (c) and (d) TEM bright-field images taken at two tilt angles illustrating (α') phase and α matrix [BAL 2011].

Figure

  Figure II.2: Heat treatment to obtain a DP microstructure: (a) schematic Fe-C diagram, (b) applied heat treatment [LAN 2011]. Although the structure of the DP steel contains ferrite and martensite, sometime a small amount of residual austenite may appear after the heat treatment. This residual microstructure could have an effect by modifying the mechanical properties. The presence of austenite phase reduces the martensite volume fraction and affects its distribution after the heat treatment. Experimental literature [SCH 2010], [KRE 2010], [NIA 2012] describes that the variation of the martensite quantity and the microstructural features of the martensite distribution affect the macroscopic behaviour of the DP steel. Referring to Tamarelli et al. [TAM 2011]observations, the ferrite is continuous for many grades up to DP780 steel, but as volume fraction of martensite exceed 50 percent (as might be found in DP980 steel or higher strengths), the ferrite may become discontinuous.

  .1 summarizes the product property requirements for various types of DP steels. It also underlines the effect of strain and bake hardening (locking dislocation by solute Carbone) [MAR 1982] [CAI 1985] [SUH 1997] [BAG 1999] [ERD 2002] [KAW 2003] [MAZ 2007] [AVR 2009].

Figure II. 4 :

 4 Figure II.4: DP steels used for safety-critical body parts (General Motors, [BUI 2011]).

Figure

  II.6: a) Springback of DP600 channel draw, b) Elastic strain recovery after the tooling is removed called Springback phenomenon [WAG 2006].

Figure II. 7

 7 Figure II.7 provides an example of the well know Forming Limits Diagrams (FLD) used to quantify formability and allows the steel and automotive makers to reduce the costs of designing tools and to shorten the time-to-market cycle. It gathers the published curves measured by different research groups (ArcelorMittal, Ramazani et al. [RAM 2012], Uthaisangsuk et al. [UTH 2008-2009].

Figure

  Figure II.8: X-Ray tomography principle [LAN 2011].

  Figure II.9(a). One specimen is pulled for a given deformation. The specimen presented in Figure II.9(b) is unloaded. The tensile device is placed on the beam line. During imaging, the deformation is stopped but maintained constant. Landron also performed a so called in-situ continuous test. It is the same procedure as the previous one but with no interruption of the deformation during the tensile test. This characterization requires a small scan time to have clear images of the reconstructed volume (as shown in Figure II.9(c)). The displacement speed of the tensile device is between 1µm/s and 5µm/s and the reconstructed volume is equal to 0.3×0.3×0.3mm³. a) b) c) Figure II.9: Experimental test setting. (a) In-situ X-ray experimental device with 1-mm notched sample [LAN 2011], (b) the specimen design, (c) Studied spatial volume for the porosity measurement at the center of the specimen, of dimensions 0.3×0.3×0.3mm³.TableII.4 presents the DP steels mechanical properties used in this study. It shows the heterogeneity of this category of AHSS in terms of stress, microstructure, and formability. DP steels with the same volume fraction of martensite for instance DPI and DPII have different material behaviour only by changing the sheet thickness or the cooling process. Indeed, DP steels can be hot or cold rolled formed (HR or CR). If hot-rolled, cooling is carefully controlled to produce the Ferritic-Martensitic structure from austenite. If continuously annealed or hot-dipped, the final structure is produced from a dual phase Ferritic-Austenitic structure that is rapidly cooled to transform some of the austenite to martensite [ARC 2012]. Krebs et al. [KRE 2010] observations on DP steel accurately shows that the hot and cold rolling processing conditions generate 'banded structures' i.e., irregular, parallel and alternating bands of ferrite and martensite, which are detrimental to mechanical properties and especially for in-use properties. Already mentioned, it is known the volume fraction of martensite has an influence on the mechanical properties. Although, DPII and DPIII have the same cooling process and sheet thickness, they have different mechanical properties. The Yield stress and the ultimate tensile stress rise by increasing significantly the volume fraction of martensite (α') appearing in the structure as seen in Figure II.10. Table II.4: Studied Dual Phase steels and their mechanical property requirement according to ArcelorMittal standard tensile test. In addition to Figure II.10 and TableII[START_REF]I.7 ~II.13 ~ performance and protect the passenger by absorbing the low and high-speed vehicle deformation over a specific distance. Beyond the safety requirement, DPI and DPIII steels give the opportunity to compare the plasticity and[END_REF] Figure II.11 shows the tensile stress-plastic strain curves of the three DP steels grades used in the study. One can easily understand the strong difference in term of tensile stress behavior in function of strain for various range of volume fraction of martensite. DPI and DPII have a good balance of strength, formability, energy absorption and durability while DPIII has a high strength and a weaker formability characteristic. Apart of this, the effect of the hard Martensitic islands is important but the role of the soft Ferrite in the final DP material is not negligible. Ferrite matrix is exceptionally ductile and absorbs strain around the Martensitic islands, enabling uniform elongation with high work hardening rate.Product Thickness [mm] YS[MPa]

Acknowledgments

II.2.2 X-ray tomography principle

The technique of X-Ray tomography started in the end of the 80's in the medical field by improving the detection of brain and breast cancer tumor. This was available by using 2D radiographies. In the middle of the 90's the X-ray tomography was adapted and found its respective place in materials science to analyze debris from plane and cars crashes for the insurance companies. Since 2000 the method has become an efficient investigation tool for various materials. The principle of X-ray tomography used by Landron and co-workers shown in Figure II.8 is composed of X-ray beam, the sample to analyze, the camera and the projections acquisition. The object to characterize rotates about a single axis while a series of 2D X-ray absorption images is recorded. Using mathematical principles of tomography, this series of images is reconstructed to produce a 3D digital image where each voxel (volume element or 3D pixel) represents the X-ray absorption at that point. Section III.2 presents how the mechanical behaviour of DP steels can be predicted with the plasticity theory. Attention will be devoted to explain the isotropic and anisotropic yield functions as well as the hardening models. Section III.3 describes the coupled damage modeling concept where mechanical behavior is affected by the damage growth due to loading. The specific case of DP steels will be investigated. Within this section, both ductile damage mechanisms and some damage models are unveiled. Preceding the conclusion, the last section is dedicated to the fracture criteria or uncoupled damage approach applied on the DP steels.

III.2 Plasticity modeling III.2.1 Yield functions

The Dual-Phase steel deforms elastically. During a monotonic loading it suddenly yields. In the plastic strain domain, the flow stress first increases due to hardening then eventually, it may soften due to damage. In numerous mechanical books The goal of this section is to define the yield function with its associative normality rule. It briefly summarizes the most commonly isotropic and anisotropic yield functions that can be used to model the plastic behaviour of the DP steels.

III.2.1.1 Yield function and the associative normality rule

The yield surface or plastic yield criterion defined in the stress space as seen Figure III.1, models the elastic limit and the beginning of the plastic flow. It is written as:

, 0

Where eqv V V is the equivalent stress and y V is the material flow stress. On one hand, when eqv V V is smaller than y V 0 p F , the deformation is purely elastic.

On the other hand when eqv

, the border is reached and the DP steel starts to plastically deform. At initial state and before hardening function takes place, the yield surface is written as following:

~III.5 ~

Where O O is the plastic multiplier.

III.2.1.2 Isotropic yield functions

The isotropic yield functions must not depend on the orientation of the load. The yield function p F is based on the deviatoric stress tensor invariants. Figure III.2 shows an overview of the most used isotropic yield criteria.

The von Mises yield criterion

The von Mises criterion [MIS 1928], known as the maximum distortion energy criterion, octahedral shear stress theory, or Maxwell-Huber-Hencky-von Mises theory, states that the plasticity is isotropic. The equation for the von Mises yield function is:

Where c V is the deviatoric stress tensor: I tr V V V


Benjamin Mathon
email: benjamin.mathon@grenoble-inp.fr

François Cayre
email: francois.cayre@gipsa-lab.grenoble-inp.fr

Member, IEEE Patrick Bas
email: patrick.bas@ec-lille.fr

Fellow, IEEE Benoît Macq
email: benoit.macq@uclouvain.be

Optimal Transport for Secure Spread-Spectrum Watermarking of Still Images

This article studies the impact of secure watermark embedding in digital images by proposing a practical implementation of secure spread-spectrum watermarking using distortion optimization. Because strong security properties (key-security and subspace-security) can be achieved using Natural Watermarking (NW) since this particular embedding lets the distribution of the host and watermarked signals unchanged, we use elements of transportation theory to minimize the global distortion (MSE). Next, we apply this new modulation, called Transportation Natural Watermarking (TNW), to design a secure watermarking scheme for grayscale images. TNW uses a multiresolution image decomposition combined with a multiplicative embedding which is taken into account at the distribution level. We show that the distortion solely relies on the variance of the wavelet subbands used during the embedding. In order to maximize a target robustness after JPEG compression, we select different combinations of subbands offering the lowest BERs for a target PSNR ranging from 35 to 55 dB and we propose an algorithm to select them. The use of transportation theory also provides an average PSNR gain of 3.6 dB on PSNR with respect to the previous embedding for a set of 2, 000 images.

Index Terms

Watermarking, Digital images, Transportation theory

I. INTRODUCTION

The development of multi-bit watermarking techniques, which consist in embedding a message (several digits, often binary symbols) in a host content such as a movie, a song, a text or a picture, has followed different paths of development since the pioneering studies. Initially, the main goal for the community was to design a watermarking scheme which can resist common signal processing operations. The embedded mark had to be detected after geometrical transformations, compressions, Gaussian noise addition, etc. The second intuitive constraint was that the presence of the watermark should remain imperceptible. These two constraints, respectively robustness and imperceptibility, are very important for a watermarking scheme but the major problem for researchers is the non-orthogonality of these constraints. In fact, in order to be robust, the mark has to be embedded with a strong power but this power implies a strong distortion of the host content. Combining robustness and imperceptibility has been a challenge for several years [START_REF] Cox | Watermarking as communications with side information[END_REF]- [START_REF] Cox | The Morgan Kaufmann Series in Multimedia Information and Systems[END_REF], which motivated the development of watermarking based on spread-spectrum modulations [START_REF] Cox | Secure spread spectrum watermarking for multimedia[END_REF], [START_REF] Hartung | Spread spectrum watermarking: Malicious attacks and counterattacks[END_REF] or quantization techniques [START_REF] Chen | Quantization index modulation methods for digital watermarking and information embedding of multimedia[END_REF]- [START_REF] Perez-Freire | Spread-spectrum vs. quantization-based data hiding: misconceptions and implications[END_REF].

A constraint receiving more and more attention in watermarking nowadays is the security constraint, defined as "the inability by unauthorized users to have access to the raw watermarking channel" by Kalker [START_REF] Kalker | Considerations on watermarking security[END_REF] and linked with the presence of an adversary [START_REF] Furon | Security analysis[END_REF]- [START_REF] Xie | Towards Robust and Secure Watermarking[END_REF]. Watermarking schemes then are required to respect the Kerckhoffs' principle [START_REF] Kerckhoffs | La cryptographie militaire[END_REF] (i.e. a secret key is shared between the encoder and the decoder and is the only secret parameter for the adversary [START_REF] Cox | Watermarking is not cryptography[END_REF], [START_REF] Bas | Practical Key Length of Watermarking Systems[END_REF]). The Kerckhoffs' principle comes from cryptography and cryptanalysis but is also widely used in watermarking [START_REF] Barni | A general framework for robust watermarking security[END_REF], [START_REF] Cayre | Kerckhoffs-based embedding security classes for woa data hiding[END_REF] or in other forms of data-hiding like steganography [START_REF] Ker | Perturbation Hiding and the Batch Steganography Problem[END_REF]- [START_REF] Moskowitz | A new paradigm hidden in steganography[END_REF]. However, watermarking security analysis shows that an adversary may estimate the secret if he owns a sufficient number of marked contents and if the embedding scheme is not secure [START_REF] Furon | Security analysis[END_REF]- [START_REF] Xie | Towards Robust and Secure Watermarking[END_REF]. According to the degree of estimation of the secret, the adversary would be able to erase, modify or copy the embedded message to another host content. There are two major differences between robustness and security. On one side, robustness attacks are deemed not intentional: signal processing operations on watermarked contents like compression can be done by the provider before the legal use by a user of this content. On the other side, security attacks come from an adversary who deliberately wants to hack the watermarking scheme, such attacks are deemed more harmful than robustness attacks because the adversary can both alter the embedded message and perform February 10, 2014 DRAFT an optimal minimization of the attack distortion [START_REF] Mathon | Comparison of secure spread-spectrum modulations applied to still image watermarking[END_REF] or even recover the original image [START_REF] Pérez-Freire | Security of lattice-based data hiding against the watermarked-only attack[END_REF].

In order to understand the ins and outs of the security game in watermarking, it is important to detail the role of the three important actors in the data hiding process:

1) the distributor is the person who marks the host contents. He is the only one who knows the secret key used for embedding and decoding. The distributor can be the author of the content, the company who sells the content or a trusted third party. His main important constraints are robustness and security,

2) the user is a person with privileges related to the content, and set by the distributor. For examples:

reading a digital book, listening to a song or watching a video. His most important constraint is the imperceptibility of the watermark,

3) the adversary will try to estimate the secret key in order to modify the watermark. Because of this modification, he can obtain more privileges than a simple user (unauthorized by the distributor). For example, in a fingerprinting scenario: he can resale the content in an illegal way, share the content on peer-to-peer networks, etc. The adversary can be a (group of) user(s) or corrupted user(s).

Robustness, imperceptibility and security are the three constraints that a distributor has to take into account for the development of a watermarking scheme. More precisely, the distributor looks for an optimization of one constraint while settling the others. Moreover, the importance of each constraint depends on the use case. The trade-off between imperceptibility, robustness and security is represented by the triangle of constraints in watermarking on Fig. 1.

In this work, we are interested in optimizing robustness and imperceptibility constraints when security is set first. The security of a watermarking scheme depends on the assumptions related to the use case under consideration. Security attacks on watermarking are mainly linked with information leakage: what are the data that the adversary can have access to?

We consider here the Watermarked Only Attack framework (WOA) [START_REF] Cayre | Kerckhoffs-based embedding security classes for woa data hiding[END_REF], which means that an adversary has access to both the source code of the scheme (without knowing the secret key) and N o ≥ 1 contents watermarked with the same, unknown key. Gathering information given by the marked contents, his goal will be to estimate the secret key, i.e. a set of secret watermark components, that were initially used during embedding. If this is not possible, the adversary can try to estimate a less informative secret parameter such as, for example, the private subspace spanned by these components in the case of spread-spectrum embedding.

It is important to point out that security in WOA depends on the distribution of contents after watermarking. We recall the three main classes of security in WOA [START_REF] Cayre | Kerckhoffs-based embedding security classes for woa data hiding[END_REF] by increasing security level:

February 10, 2014 DRAFT Figure 1. The triangle of constraints in watermarking where the distributor looks for an optimization of one constraint while settling the two others. The importance of each constraint depends on the use case scenario and the distributor has to take into account the trade-off between the three constraints.

• Insecurity: the distribution of marked contents is unique and specific to the chosen secret key. More precisely, two different keys generate two different distributions of marked contents. An adversary can then estimate the secret components [START_REF] Mathon | Practical performance analysis of secure modulations for woa spread-spectrum based image watermarking[END_REF].

• Key-security: for a subset of keys, the distribution of the marked contents will be the same. With a sufficient number of observations, the adversary will be able to estimate a subset of the watermarking space (a subspace in the case of spread-spectrum watermarking).

• Subspace-security: the distribution of the marked contents is the same for each possible key. The adversary cannot estimate a set smaller than the original set of possible watermark signals.

Spread-spectrum modulations such as Natural Watermarking (NW) [START_REF] Bas | Natural Watermarking: A Secure Spread Spectrum Technique for WOA, ser[END_REF] and Circular Watermarking (CW) [START_REF]Achieving subspace or key security for woa using natural or circular watermarking[END_REF] have been designed to deter a correct estimation of the secret components: NW prevents the whole subspace estimation and is consequently subspace-secure, while CW prevents the estimation of the watermark components and therefore qualifies for key-security.

In this article, we propose a practical implementation of NW for still images. We look for an optimization of the distortion constraint that reduces the distortion due to the assignment of each host content to a marked content. This problem finds a solution using transportation theory [START_REF] Villani | Optimal transport: old and new[END_REF], a research domain that was successfully used for the development of watermarking techniques as in [START_REF] Mathon | Optimization of natural watermarking using transportation theory[END_REF] (Transportation NW modulation), [START_REF] Bas | Soft-SCS: improving the security and robustness of the Scalar-Costa-Scheme by optimal distribution matching[END_REF] (Soft-scalar Costa Scheme) and [START_REF] Cao | Controllable secure watermarking technique for tradeoff between robustness and security[END_REF] (Controllable Secure Watermarking). We implement the Transportation NW modulation (TNW) on multiple combinations of wavelet subbands of grayscale images and we add a perceptual weighting in order to increase the imperceptibility of the February 10, 2014 DRAFT whole system. We show that our optimization of imperceptibility can be made while keeping the desired security level. This work is a practical extension of [START_REF] Mathon | Optimization of natural watermarking using transportation theory[END_REF], where the transportation problem was presented and applied only to synthetic Gaussian signals, and we focus here on the design of secure embedding on still images.

In Sec. II, we illustrate the links between subspace-secure watermarking and the distributions of marked contents using two embedding techniques. We present elements of transportation theory which optimize assignments between host and marked distributions in a theoretical way, and we further use these results

to design TNW. In Sec. III, we tackle implementation issues and propose an experimental watermarking scheme in the wavelet domain with psycho-visual masking. Since the robustness and the distortion of the proposed scheme are intimately linked with the set of the selected subbands, we present the combinations which provide the best robustness for a given distortion, computed from 2, 000 images.

II. DISTRIBUTION MATCHING WATERMARKING

We first list the notations and conventions used in this article. Functions are denoted in roman fonts, sets in calligraphy fonts, vectors and matrices in bold fonts and variables in italic fonts. Vectors are written in small letters and matrices in capital ones. The content of a vector x with length n is denoted by (x(0) . . . x(n -1)), the random variable associated to i.i.d components is denoted X. P δ (resp. f δ )

is the cumulative distribution function (resp. probability density function) of a distribution δ. |.| is the absolute value on R and . is the Euclidean norm on R N with N ≥ 1. .|. is the usual inner product.

A. Settling the subspace-security constraint

In this work, we are focusing on the subspace-security class of the WOA framework. Following Kerckhoffs' principle, an adversary is deemed to have access to all public data of the watermarking scheme: transformation domain where the embedding takes place (e.g. DCT, DWT) and can therefore access the feature vectors of size N v of the contents. For the definition of the secret key, we use the same formalism as in [START_REF] Cayre | Kerckhoffs-based embedding security classes for woa data hiding[END_REF], [START_REF] Costa | Writing on dirty paper[END_REF]: a secret key K s generates a set of N h secret components, C = {c i } i∈Nh . A message m ∈ {0, 1} Nc is associated to one or more secret components, meaning that N h ≥ 2 Nc .

We denote:

• X = (x 0 . . . • S: the set of embedding and decoding keys.

To efficiently modify the embedded message, the first step for an adversary is to estimate the secret components of K s . According to Kerckhoffs' principle, he has access to:

• p(X): the distribution of host feature vectors which can be explicitly known (analytic formula) or modeled (the adversary can generate his own signals with contents he may have collected elsewhere),

• p(Y|K i ): the conditional distribution of marked signals given a key K i ∈ S, K i being generated by the adversary,

• p (Y Ks ): the distribution of signals marked by the distributor with the secret key K s , the adversary has access to the marked contents without the knowledge of secret key K s .

A watermarking scheme qualifies for subspace-security if the distribution of marked signals is the same for each possible key. More precisely, the pdf p(Y|K i ) does not depend on the key K i . This assumption is equivalent to state that Y and K i are independent. Formally [START_REF] Cayre | Kerckhoffs-based embedding security classes for woa data hiding[END_REF]:

∀K i ∈ S, p (Y|K i ) = p (Y Ks ) . (1)
In this case, an adversary cannot estimate K s because he sees no differences by trying several secret keys in S. One way for a distributor to achieve subspace-security in WOA framework is to ensure that the distribution of the watermarked signals is identical to the distribution of the original host signals.

Note that this condition implies that the watermarking scheme also belongs to the key-security class.

It means that there exists a subset of keys S c ⊂ S where distributions of marked signals are the same.

An adversary can estimate this subset of keys which contains the secret key K s :

∀K i ∈ S c , p (Y|K i ) = p (Y Ks ) . (2)
Finally, a watermarking scheme belongs to the insecurity class if for each secret key, the distribution of marked contents is unique, formally:

∃! K i ∈ S, p (Y|K i ) = p (Y Ks ) , (3)
and:

∀K 1 = K 2 ∈ S, p (Y|K 1 ) = p (Y|K 2 ) . (4)
February 10, 2014 DRAFT

We now review some previous relevant works on secure WOA watermarking.

1) χ 2 Watermarking: In [START_REF] Mathon | Distortion Optimization of Model-Based Secure Embedding Schemes for Data-Hiding[END_REF], we have proposed an illustrative subspace-secure scheme called χ 2

Watermarking (χ 2 W). This technique modifies the square Euclidean norm of a host signal to embed a message, and the secret components are then partitioned along the positive real axis R + according to a secret key. In fact, if x ∈ R Nv with x ∼ N (0, 1),

x 2 ∼ χ 2 (N v ) (χ 2 law of degrees N v ).
The secret components are then set in a partition of [0, +∞[. To embed a secret message m in a host vector x, we randomly choose a square norm N(m) in the corresponding real interval and we compute:

y = N(m) x 2 x. (5)
Watermarking is then done by selecting a square norm in the secret component corresponding to the message. By keeping the distribution of the square norms after watermarking, this technique is subspacesecure. Fig. 2 shows the distribution of the norms of 2, 000 watermarked signals y j 2 with messages of N c = 2 bits.

However, a severe drawback of χ 2 W is its weak robustness against white Gaussian noise addition, contrary to other methods such as spread-spectrum modulations.

2) Spread-spectrum Watermarking: Spread-spectrum techniques are good candidates for robust watermarking, which can also be extended to secure versions. We consider a message m ∈ {0, 1} Nc to be embedded into a host Gaussian distributed content x j ∈ R Nv . The message is coded using N c carriers u i ∈ R Nv i∈Nc which are the secret components. These carriers are generated with a pseudo random number generator (PRNG) initialized with a seed K s ∈ N. They come as zero-mean Gaussian vectors obtained with the PRNG and are further orthogonalized to provide a basis of the private subspace (Sec. II-A4 provides more details about the generation of the carriers). In the WOA framework, security attacks are connected with the estimation of the carriers u i . It is not necessary to go back to the PRNG key K s to perform a security attack. For the rest of this article, we make no difference between K s and the carriers u i .

For each host signal x j , the embedding function e Ks (depending on the secret key K s ) creates the marked signal y j (in order to hide the message m j ) following:

∀j ∈ [N o ], y j = e Ks (x j , m j ) = x j + w j = x j + Nc-1 i=0 s(m j (i), x j )u i , (6)
where w j is the watermark signal and s :

{0, 1} × R Nv → R is a modulation. Decoding is performed using correlations ỹj ∈ R Nc : ∀i ∈ [N c ], ỹj (i) = 1 N v y j |u i . (7)
We obtain:

February 10, 2014 DRAFT mj (i) =      0 if ỹj (i) > 0, 1 if ỹj (i) ≤ 0, (8)
where mj (i) is the i-th decoded bit.

The classical Spread-Spectrum modulation (SS) is defined by:

s SS (m j (i), x j ) = α(-1) mj(i) , ( 9
)
where α is a parameter driving the strength of the embedding, i.e. the distortion.

This modulation is known to belong to the insecure watermarking class [START_REF] Cayre | Kerckhoffs-based embedding security classes for woa data hiding[END_REF], as secret carriers can be precisely estimated using first Principal Component Analysis (PCA) to estimate the secret subspace then Independent Component Analysis (ICA) [START_REF] Mathon | Comparison of secure spread-spectrum modulations applied to still image watermarking[END_REF], [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF] to estimate the secret components (see II-A3).

Subspace-secure spread-spectrum watermarking has been proposed in [START_REF] Bas | Natural Watermarking: A Secure Spread Spectrum Technique for WOA, ser[END_REF] with the Natural Watermarking (NW), this modulation enables to keep the same Gaussian distribution of correlations between host signals and secret carriers.

The s N W modulation is given by:

s N W (m j (i), x j ) = (-1) mj(i) sign(x j (i)) -1 xj (i). (10)
This technique uses symmetries of host correlations between the N c secret carriers to embed the N c bits with classical spread-spectrum decoding rules by keeping the relation |ỹ j (i)| = |x j (i)|.

3) Practical security analysis: Considering that an adversary has access to N o spread-spectrum watermarked contents Y Ks , his goal is to obtain as much information as possible about the secret components U = (u 0 . . . u Nc-1 ). Principal Component Analysis allows for an adversary to estimate the N cdimensional private subspace spanned by the secret components if the embedding alters the covariance matrix of the contents. This technique aims at finding the optimal linear transformation corresponding to the subspace yielding the largest variance. For classical SS modulation (Eq. ( 9)), message embedding increases the variance of the signal in the direction of the carriers. Note that PCA deals with subspacesecurity: contrary to classical SS technique, this technique cannot be correctly applied to contents watermarked with theoretical NW modulation. In Sec. III-D, we analyse the security of NW modulation used in a practical implementation (for still images) using PCA.

February 10, 2014 DRAFT

To measure the accuracy of the estimation of the private subspace, we use the chordal distance [START_REF] Pérez-Freire | Spread-spectrum watermarking security[END_REF], [START_REF] Conway | Packing lines, planes, etc.: Packings in grassmannian spaces[END_REF], which provides a distance between two subspaces. If Û = (û 0 . . . ûNc-1 ) denotes the estimated carriers, the chordal distance between U and Û is defined by:

d c = 1 √ N c Nc-1 i=0 sin 2 (θ i ) 1/2 , ( 11
)
where the θ 0 . . . θ Nc-1 denote the principal angles between U and Û. The chordal distance d c equals 0 when the subspaces spanned by the both matrices are the same and 1 when the matrices are orthogonal.

4) On the generation of the secret carriers: While there is no explicit agreement in the watermarking community on how to generate secret carriers, we describe here how we precisely generate ours:

• The Mersenne Twister MT19937 [START_REF] Matsumoto | Mersenne twister: A 623-dimensionally equidistributed uniform pseudorandom number generator[END_REF] PRNG is seeded with the secret key (this PRNG has no known bias, generates high-quality output for simulations and accepts seeds as long as 768 bits). We obtain uniformly distributed vectors;

• Going from a uniform to a Gaussian distribution is done with the Ziggurat method [START_REF] Marsaglia | The ziggurat method for generating random variables[END_REF]. This method is more reliable and quicker than the usual Box-Muller transform [START_REF] Box | A note on the generation of random normal deviates[END_REF];

• In order to remove inter-symbol interference, we explicitly orthogonalize the vectors with the help of the Gram-Schmidt procedure (we have verified that this procedure has a very limited impact, if any, on the gaussianity of the vectors). Further, we also explicitly remove the mean and normalize the output to obtain N (0, 1)-distributed vectors. These vectors are now the secret carriers we use.

In itself, the MT19937 is not a cryptographically secure PRNG. It is vulnerable to a Berlekamp-Massey attack [START_REF] Berlekamp | Nonbinary BCH decoding[END_REF], [START_REF] Massey | Shift-register synthesis and bch decoding[END_REF]. However, performing such an attack in the context of digital watermarking has not been done before and would prove extremely difficult because:

• An adversary would need to recover the exact, binary output of the MT19937 to perform the attack;

• He could only estimate the normalized, zero-mean output of the Ziggurat algorithm. To the best of our knowledge, the Ziggurat algorithm is not invertible and even the quality of the adversary's estimation may be insufficient in practice (when such an estimation is possible, i.e. when the embedding technique at hand is insecure). Also, the adversary would need to guess the original mean and variance of the vectors after applying the Ziggurat algorithm, all of which is lost information after the final step.

However, in case such an attack would prove feasible in practice, which we doubt, it would be sufficient to change for a cryptographically secure PRNG based, e.g., on a standard cryptosystem. Thus, carriers estimation is the only attack one should consider in the context of digital watermarking security.

February 10, 2014 DRAFT

B. Optimization of imperceptibility

For the two subspace-secure techniques χ 2 W and NW we have presented, in order to keep the distribution of the correlations after watermarking, it is not possible to change the distortion after embedding.

Security of watermarking schemes in WOA framework can be completely defined by the distribution of signals after embedding (distribution of correlations for spread-spectrum techniques, distribution of Euclidean norms for χ 2 W), and this property therefore constraints the distortion.

In this Section, watermarking in WOA is presented from a new perspective: we first generate a marked distribution in the subspace spanned by the secret key (a bijection between the security level one wants to achieve and the distribution). Next, we create marked signals by retro-projecting into the N v -dimensional space.

From the distortion point of view, the subspace-secure techniques presented above are not optimal.

In the WOA framework, we want to watermark N o host contents with random messages. For each host content (a point in the host distribution in the private subspace), we associate a marked content in the targeted distribution. The question now becomes: how to associate each point of a host distribution to each point of a marked distribution in a decoding region by minimizing the distances between these points (the distance being proportional to the distortion induced by the watermark)?

1) The Hungarian optimization algorithm: In [START_REF] Mathon | Distortion Optimization of Model-Based Secure Embedding Schemes for Data-Hiding[END_REF], we have proposed a way to match a host distribution in the private subspace with a marked distribution, while minimizing the global distance using the Hungarian algorithm [START_REF] Kuhn | The hungarian method for the assignment problem[END_REF]. This method is discrete and it is based on a precomputed matching on N m host correlations (N m ≥ N o ) vectors X and N m marked correlations vectors Ỹ. The Hungarian algorithm is able to find the bijection H which minimizes the Euclidean distance (weight) on average between X and Ỹ. To watermark a signal x, we find the nearest neighbor of x in X and we obtain ỹ based on the optimal matching given by the Hungarian algorithm. However, this method has three disadvantages. First, the complexity is O(N m 3 ) and can be a real computational problem when N m and N o are important.

Secondly, the performances (gain of distortion) decrease when the dimension of the private subspace (the number of secret components for SS techniques) increases, this is due to the difficulty to find Nearest Neighbors in high dimensional spaces. Finally, this method is bounded by the number of observations we use to compute the matching and it becomes largely sub-optimal for points lying beyond the tail of the empirical distribution (see Fig. 5). Fig. 3 shows an illustration of two subsets X and Ỹ (N m = 12), as well as the optimal matching H found by the Hungarian algorithm (weights are not depicted).

To avoid these problems, we present in the next subsection a method based on the transportation theory which minimizes the global distortion and its application on subspace-secure spread-spectrum watermarking.

2) Transportation theory: This theory has been first considered by mathematicians Monge [START_REF] Monge | Mémoire sur la théorie des déblais et des remblais[END_REF] and Kantorovitch [START_REF] Kantorovitch | On the translocation of masses[END_REF], [START_REF] Kantorovich | On a problem of monge[END_REF]. This domain consists in finding a bijection T * (a transport map) in order to move an initial set of points X ⊂ R N with distribution µ to a final set Ỹ ⊂ R N with distribution ν given a cost function c : X × Ỹ → R + . T * is an optimal transport map if it minimizes the total cost as illustrated in Fig. 4. Formally, the problem is to find:

T * = arg min T ˆRN c(x, T(x))µ(x)dx ν = µ • T -1 . (12)
In [START_REF] Villani | Optimal transport: old and new[END_REF], [START_REF] Rachev | Mass Transportation Problems[END_REF], if h, defined by h(x -ỹ) = c(x, ỹ), is a convex function, a one-dimensional solution (N = 1) is given by:

T * = P -1 ν • P µ . (13)
For N ≥ 1, we can apply the Knott-Smith criterion [START_REF] Knott | On the optimal mapping of distributions[END_REF]: if c(x, ỹ) = x-ỹ 2 , the sufficient conditions implying that T * is a minimizer of Eq. ( 12) are: i) T * ( X ) ∼ ν (image of X by T * is distributed according to ν),

ii) the Jacobian matrix J T * of T * is symmetric and positive semi-definite. One way to achieve subspace-security in spread-spectrum watermarking is to keep the distribution of the correlations of the contents over the N c secret carriers. One can see the links between transportation theory and secure watermarking. Given a host distribution, it is possible to generate a marked distribution while minimizing the distortion cost. In our problem, we consider the minimization of the quadratic distance, in order to maximize the PSNR for instance, hence we set c(x, ỹ) = xỹ 2 .

3) Optimal embedding: We denote with X = {x j } j∈No ⊂ R Nc the set of host correlations and Ỹ = {ỹ j } j∈No ⊂ R Nc is the set of watermarked correlations. Further, we assume that host contents x are Gaussian distributed. This assumption is plausible due to the central limit theorem and to the fact that we work on host components projected on random vectors. We have therefore:

∀i ∈ [N c ], x(i) ∼ N 0, σ 2 x /N v = µ, (14)
with:

P µ (t) = 1 2 1 + erf t √ N v σ x √ 2 . ( 15
)
Since spread-spectrum correlations are symmetric in the private subspace, we first consider the embedding of a N c -bit constant message for each host content, for example m = 0 = (0 . . . 0), whose decoding February 10, 2014 DRAFT region is represented by the region R + Nc in the private subspace. To guarantee subspace-security, we need:

∀i ∈ [N c ], ỹ(i) ∼ N + (0, σ 2 x /N v ) = ν, (16)
where N + denotes a Gaussian distribution truncated in R + .

If δ ∼ N 0, σ 2 x /N v , ν ∼ N + (0, σ 2 x /N v ) is given by its cdf:

P ν (t) =      0 if t < 0, 2P δ (t) -1 if t ≥ 0. (17)
The quantile function of P ν is given by:

P -1 ν (t) = P -1 δ (t/2 + 1/2) = σ x √ 2 √ N v erf -1 (t). (18)
The strategy we adopt here is the following: we apply an optimal transport map for each dimension of the N c -dimensional subspace using Eq. ( 13), Eq. ( 15) and Eq. ( 18). We consequently have the transport map T 0 : R Nc → R Nc :

ỹ = T 0 (x) =      P -1 ν • P µ (x(0)) . . . P -1 ν • P µ (x(N c -1))      . ( 19
)
Proposition 1. T 0 is an optimal transport map.

Proof:

This proposition is intuitively straightforward due to the separability of the problem along the N c orthogonal secret components. However, in order to prove it rigorously, we have to use the Knott-Smith criterion (Sec. II-B2) to prove that T 0 is an optimal transport map. Assumption i) is verified because ỹ is a Gaussian vector (its components are Gaussian i.i.d.). In order to verify the assumption ii), we introduce the Jacobian matrix of T 0 :

February 10, 2014 DRAFT

J T0 =         ∂(P -1 ν •Pµ(x(0))) ∂ x(0) 0 • • • 0 0 . . . . . . . . . . . . . . . . . . 0 0 • • • 0 ∂(P -1 ν •Pµ(x(Nc-1))) ∂ x(Nc-1)         . ( 20
)
We now prove that this matrix is positive semi-definite. By construction of T 0 , J T0 is symmetric.

Then J T0 is positive semi-definite iff eigenvalues of J T0 are positives or null. We now have to prove that P -1 ν • P µ (t) ≥ 0. We have:

P -1 ν • P µ (t) = P -1 ν (P µ (t)) f µ (t). ( 21
)
f µ is positive (pdf) and P -1 ν is a non-decreasing function, so P -1 ν (t) ≥ 0. Finally P -1 ν • P µ (t) ≥ 0 and J T0 is positive semi-definite: T 0 is an optimal transport map.

Because of the symmetry property of decoding with correlations in spread-spectrum watermarking, in order to embed messages which differ from 0, central symmetries have to be performed on components of x and ỹ before and after embedding to correctly hide the message m. For any m ∈ F Nc 2 , the optimal transport map T m is given by:

T m (x) = R -1 (T 0 (Rx)), (22)
with:

R ∈ M Nc,Nc (R), R(i, j) =    0 if i = j, (-1) m(i) otherwise. ( 23
)
This optimal transport map allows us to design a new NW modulation called Transportation Natural Watermarking (TNW) where the modulation is given by:

s T N W (m(i), x) = T m (x)(i) -x(i). ( 24
)
We can also compute the distortion gain offered by TNW over NW, e.g. the difference between the two Watermark to Content Ratios W CR = 10 log 10 (σ 2 w /σ 2 x ). Without loss of generality, we compute

E [W CR] for s = s N W and s = s T N W considering m(i) = 0.
February 10, 2014 DRAFT For s = s N W :

E s(m(i), x) 2 = E ((x(i)/|x(i)| -1) x(i)) 2 = 1/2×E 4x(i) 2 = 2/N v 2 E x|u i 2 = 2/N v 2 E     Nv-1 j=0 x(j)u i (j)   2   = 2/N v 2 × Nv-1 j=0 E x(j) 2 × E u i (j) 2 = 2σ 2 x /N v . ( 25
)
We obtain then:

E [W CR N W ] = 10 log 10 σ 2 w σ 2 x = 10 log 10 2N c N v . (26)
For s = s T N W :

E s(m(i), x) 2 = ˆ+∞ -∞ P -1 ν • P µ (t) -t 2 f µ (t)dt. (27)
A = E s(m(i), x) 2 has to be computed using numerical integration. We obtain:

E [W CR T N W ] = 10 log 10 σ 2 w σ 2 x = 10 log 10 A × N c σ 2 x . (28)
The WCR gain using TNW modulation instead of NW modulation is given by:

Gain [dB] = E [W CR N W ] -E [W CR T N W ] = 10 log 10 2σ 2 x AN v . ( 29
)
After performing numerical integration, we notice that

E [W CR T N W ] is constant w.r.t σ 2
x and that the gain is constant w.r.t N c . This gain is approximately ≈ 3.77 dB with N v = 256. Fig. 5 shows the embedding functions e (Eq. ( 6)) for NW, TNW and Hungarian HNW [START_REF] Mathon | Distortion Optimization of Model-Based Secure Embedding Schemes for Data-Hiding[END_REF] modulations functions of scalar host signal x ∼ N (0, 1), where only one bit is embedded. We can see that, unlike classical modulation, improved methods are not linear. Moreover, because HNW is based on nearest neighbors on precomputed Gaussian mapping (N m = 1, 000 in this example), this method deals with edge effects (the min and max values) when x is increasing. Fig. 6 shows the distortion (measured by the WCR) caused by watermark embedding for NW, TNW and HNW modulations and its evolution w.r.t. the number of inserted bits N c . As expected, the distortion gain between NW and TNW is constant w.r.t. the length of the embedded messages, these practical computations also match the theoretical derivations in ( 26) and ( 28) by verifying a constant gap. Additionally, we can notice that the distortion gap between NW and HNW decreases when the dimension of the private subspace (N c ) increases, this is due to the fact that the nearest neighbor search involves larger Euclidean distances and consequently incurs a loss of efficiency in high dimensional spaces. and HNW decreases when the dimension of the private subspace (Nc) increases because of larger Euclidean distances in large dimensions. On the contrary, TNW does not suffer from dimension issues and produces better results than HNW.

III. APPLICATION ON STILL IMAGES

A. From theoretical to practical secure watermarking Most watermarking schemes assume a Gaussian distribution of host signals. However, this model does not fit usual distributions of image components. For example, DCT components are often modeled by a Laplace distribution and DWT components by a generalized Gaussian distribution. In order to be close to the Gaussian assumption, we adopt the following strategy: we project the host feature vectors onto pseudorandom carriers generated with the same distribution. Following the Central Limit Theorem, the projected February 10, 2014 DRAFT signals (marginal distributions) are asymptotically Gaussian. By working on Gaussian distribution and applying NW, we can guarantee that we are subspace-secure in the projected subspace. If one is able to estimate the projected subspace, using denoising techniques for example [START_REF] Wiener | Extrapolation, interpolation, and smoothing of stationary time series: with engineering applications[END_REF], then the key-security property still holds because of the security of the Gaussian distribution created by the NW embedding.

B. Experimental proposed scheme

Fig. 7 presents our experimental watermarking scheme on grayscale images. After a 5-level 9/7

Daubechies DWT transform [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF] on the host image, we select N t components on 9 subbands in low/mean frequencies as presented in Fig. 8 to obtain a PSNR between original and marked images between 35 and 55 dB on average. Note that this DWT transform is not an orthogonal but a biorthogonal transform [START_REF] Usevitch | Optimal bit allocation for biorthogonal wavelet coding[END_REF],

however the 9/7 filter set deviates by only a few percent from orthogonal filter weighting [START_REF]A tutorial on modern lossy wavelet image compression : Foundations of jpeg 2000[END_REF]. The extracted feature vector with size N t is denoted as x t . In order to respect a normal distribution on host contents to be able to apply TNW, we use a projection of the feature vector on

N v carriers {a i } i∈Nv with size N t : ∀i ∈ [N v ], x(i) = Nt-1 j=0 x t (j)a i (j). (30)
Vectors a i are pseudo-randomly generated according to an uniform distribution: ∀i ∈ N v , ∀j ∈ N t , a i (j) ∼ U(-3/N c , 3/N c ), quasi-orthogonal and

E[||a|| 2 ] = 1. For each i ∈ N v , x(i) is
asymptotically Gaussian distributed due to the CLT. Note that the independence between each component is not provided. However, in this work, we assume this condition, partially justified by the important length of the considered feature vectors.

The watermark signal w (Eq. ( 6)) is constructed by spread-spectrum TNW modulation for a message m and N c secret carriers u i . To preserve the host distribution on the projected space and apply the TWN modulation [START_REF] Bas | Practical Key Length of Watermarking Systems[END_REF], we first estimate values of σ 2

x for each selection of subbands on the whole BOWS2-IMAGES database [START_REF] Bas | Break Our Watermarking System 2nd edition[END_REF] which contains 10, 000 images.

The watermark is then computed in the wavelet domain using the following inverse projection:

∀i ∈ [N t ], w t (i) = Nv-1 j=0 w(j)a j (i). (31)
In order to improve imperceptibility, we propose to add a psychovisual masking proposed by Piva [START_REF] Piva | DCT-Based Watermark Recovering Without Resorting to the Uncorrupted Original Image[END_REF] once the watermark is computed in the wavelet domain (Eq. ( 31)). Multiplicative embedding is consequently used to compute the final watermarked vector y t :

February 10, 2014 DRAFT Figure 7. Experimental image watermarking scheme.

y t = x t + w t with ∀i ∈ [N t ], w t (i) = 1 E[|X t |] |x t (i)|w t (i). (32)
The normalization factor (E[|X t |]) -1 guarantees that the distributions related to the projections w t |a j and w t |a j remain identical after the psychovisual weighting (E[

W t |A j ] = E[ W t |A j ]
). Note that such a strategy does not change the optimality of the embedding w.r.t. the quadratic distance since in the secret subspace it is equivalent to the identity application.

February 10, 2014 DRAFT

C. Subbands selection for distortion and decoding issues

Natural Watermarking (NW) and our improved Transportation Natural Watermarking (TNW) do not allow to set the embedding distortion in the N v -dimensional subspace. However, the strength of the watermark can be tuned by selecting appropriate combinations of subbands in low-frequencies components.

Fig. 8 depicts the 9 different components that can be selected to perform the embedding and we look for the selection which offers both the desired distortion and best robustness. We experiment here our scheme on the 2 9= 512 possible combinations of subbands on 2, 000 512 × 512 natural images from the BOWS2-IMAGES database [START_REF] Bas | Break Our Watermarking System 2nd edition[END_REF]. We hide here N c = 16 bits on each image and set the length of projected signals with N v = 256. To test the robustness, we apply a JPEG compression with quality factor Q F = 30 after watermarking and before decoding.

Fig. 9 shows the histogram of PSNR, on average, for the different combinations. The average PSNR is 44.05 dB with a deviation of 2.64 dB and the PSNRs range between 35 dB and 55 dB. Note that the relation between PSNR and WCR is given by [START_REF] Mathon | Comparison of secure spread-spectrum modulations applied to still image watermarking[END_REF]:

P SN R = 10 log 10 255 2 σ 2 x × 512 2 N v × E[|X t |] 2 E[X 2 t ] -W CR. (33)
Fig. 10 (BER, PSNR) which offer the lowest BER are depicted by square marks.

Several remarks can be outlined from these results:

1) The PSNR gain obtained using TNW instead of NW modulation is 3.6 dB on average with a standard deviation of 0.21 dB on the 512 combinations. This gain matches our theoretical computation using Eq. ( 29) and Eq. (33).

2) The lowest PSNRs correspond to low-frequencies subbands and the largest ones to mean-frequencies subbands.

3) The couples (BER, PSNR) with the lowest BERs correspond to combinations of contiguous or nearly contiguous sets of subbands. In order to increase the robustness for a given distortion budget, the best strategy is to pick the maximum number of low frequency subbands (in order to increase the spreading factor). using the same [1 2] combination for NW and TNW modulations. Fig. 13 presents a zoom of upper areas of these images. We show here the consequences of using TNW instead of NW modulation and multiplicative embedding (Eq. ( 32)) instead of constant embedding (Eq. ( 31)) from the imperceptibility point of view. Fig. 12 shows the robustness evaluation of our experimental watermarking scheme against JPEG compression: Bit Error Rate functions of JPEG Quality Factor (Q F ). We present here the four combinations which give the smallest BER after JPEG compression with Q F = 30.

D. Security assessments

We assess the security analysis of the proposed scheme using two means:

1) Distributions: we perform visual inspection in the subspace spanned by the secret components in order to verify that the distribution after embedding is not modified. 30)) and multiplicative embedding (Eq. ( 32)).

2) PCA: we try to estimate the secret subspace using PCA and measuring the distance between the estimated space and the true one using the chordal distance, then we compare our results with traditional, insecure SS. We apply here the method introduced in Sec. II-A3 It consists in estimating, using Principal Component Analysis (PCA) on the watermarked contents in R Nv , the subspace spanned by the secret components. Due to the embedding process, the variance of this subspace increases for insecure For TNW, we obtain a chordal distance of 0.98 using the N c first eigenvectors, which means that PCA is not able to correctly estimate the secret subspace. For SS, the chordal distance equals 0.08 when we focus on the N c first eigenvectors which means that in this case the estimated and secret subspaces are very close. We can then conclude this analysis by stating that the proposed implementation is immune to Distributions keep Gaussian after projection on quasi-orthogonal uniform carriers (Eq. ( 30)) and multiplicative embedding (Eq. ( 32)).

security attacks based on subspace estimation: this is due to the subspace-security property of the TNW embedding.

IV. CONCLUSION

This paper shows how to consider image watermarking by taking into account the three fundamental constraints: the security constraint is set to a given class (key or subspace-security), the distortion is minimized, and the robustness is maximized.

We have proposed to minimize the embedding distortion using both optimal transportation theory to guarantee the security class, and multiplicative embedding to minimize the visual impact. In the context of secure embedding, the robustness solely relies on the variance of the components of a given image and it is maximized by selecting the appropriate subbands. These two optimizations enable on one hand to obtain a PSNR gain of 3.6 dB on average, and on the other hand to provide the sets of configurations which maximize the robustness for a given distortion. Such a methodology can be also applied on other classes of secure watermarking schemes, by using optimal transport whenever it is possible [START_REF] Bas | Soft-SCS: improving the security and robustness of the Scalar-Costa-Scheme by optimal distribution matching[END_REF], [START_REF] Cao | Controllable secure watermarking technique for tradeoff between robustness and security[END_REF] or finding approximations of the optimal mappings for more complex distributions such as circular watermarking [START_REF]Achieving subspace or key security for woa using natural or circular watermarking[END_REF].

x

  No-1 ): the set of N o host original signals of size N v (extracted from N o host contents), • M = (m 0 . . . m No-1 ): the set of N o binary messages of size N c (in the WOA framework, we assume that the digits of the messages are chosen at random), February 10, 2014 DRAFT • Y Ks = (y 0 . . . y No-1 ): the set of N o signals marked with the secret key K s . Each signal y j hides the message m j on x j with j ∈ [N o ],

February

Figure 2 .

 2 Figure 2. Distribution of square Euclidean norms of watermarked signals. We keep the χ 2 (Nv) distribution of host norms. Parameters: No = 2, 000, Nc = 2, Nv = 55.

FebruaryFigure 3 .

 3 Figure 3. Example of triplet ( X , Ỹ, H) with Nm = 12. H is the bijection obtained by the Hungarian method which minimizes the sum of weights of edges with complexity O(Nm 3 ).

February

Figure 4 .

 4 Figure 4. Illustration of the transportation problem: we look for an optimal transport map T * : X → Ỹ with X ∼ µ and Ỹ ∼ ν which minimizes the cost c : X × Ỹ → R + .

FebruaryFigure 5 .

 5 Figure 5. Embedding functions e with Nv = 1 and m(0) = 0. The Hungarian algorithm is run with Nm = 1, 000. We can see that the mapping provided by the Hungarian method is inaccurate for values above the maximum value of the empirical distribution.

FebruaryFigure 6 .

 6 Figure 6. Distortion (WCR) w.r.t. the number of bits Nc for NW, HNW and TNW modulations, Nv = 512, σ 2 x = 1. The Hungarian algorithm is computed with Nm = 10, 000 Gaussian signals and the distortion is computed on average on No = 2, 000 Gaussian signals for the three modulations. As mentioned in Sec. II-B1, the distortion difference between NW

Figure 8 .

 8 Figure 8. Component selection: the watermark is hidden on each possible selection using 9 subbands in low-frequencies components after a 9/7 Daubechies DWT with a 5-level decomposition in order to obtain a PSNR between 35 and 55 dB.

4 )FebruaryFrequencyFigure 9 .

 49 Fig.11presents one of the 2, 000 images of 512×512 pixels watermarked with our experimental scheme

FebruaryFigure 10 .

 10 Figure 10. Average embedding PSNR w.r.t. the average BER after JPEG compression (QF = 30) on 2, 000 still images for 2 9 = 512 possible combinations. Square plots show the best combinations (low BER) for PSNR between 35 and 55 dB and cross plots are selected according an iterative procedure on subbands as presented in Tab. I.

FebruaryFigure 11 .Figure 12 .

 1112 Figure 11. Host image (a), Watermarked image using the NW modulation with constant embedding (CE) (b), the TNW modulation with CE (c) and the TNW modulation with multiplicative embedding (ME) (d). Parameters : combination = [1 2], Nt = 512, Nv = 256, Nc = 16.

FebruaryFigure 13 .Figure 14 .

 1314 Figure 13. Zoom of upper areas of images presented in Fig. 11. Host image (a), Watermarked image using NW with constant embedding (CE) (b), TNW with CE (c) and TNW with multiplicative embedding (ME) (d). Parameters : combination = [1 2], Nt = 512, Nv = 256, Nc = 16.

  presents the relation between the PSNR computed over the original and watermarked images, and the Bit Error Rate after JPEG compression with Q F = 30 for the 512 combinations. The couples COUPLES (BER, PSNR) (IN CROSS PLOTS IN FIG. 10) SELECTED ACCORDING TO AN ITERATIVE PROCEDURE WHICH CONSISTS, IN A CONTINUOUS WAY, IN SELECTING HLi AND LHi SUBBANDS (AND HHi TO IMPROVE ROBUSTNESS) AND REMOVING LOW ONES UNTIL THE DESIRED DISTORTION IS ACHIEVED.

	PSNR [dB] BER	Combination
	35.76	6.509375e-02 [1 2]
	36.56	5.478125e-02 [1 2 3]
	39.53	4.034375e-02 [1 2 3 4 5]
	39.40	4.121875e-02 [1 2 3 4 5 6]
	43.07	5.334375e-02 [4 5]
	44.05	6.268750e-02 [4 5 6]
	46.97	8.412500e-02 [4 5 6 7 8]
	47.56	8.821875e-02 [4 5 6 7 8 9]
	50.03	1.302812e-01 [7 8]
	51.01	1.209375e-01 [7 8 9]
		Table I
	February 10, 2014		DRAFT

February 10, 2014 DRAFT

This work was supported in part by the National French project ANR-10-CORD-019 Estampille, the Belgian project BCRYPT IAP (PAI) phase-VI and by the installation grant of Prof. S. Voloshynovskiy, University of Geneva, Switzerland.