Angelo Favini

Rabah Labbas

Keddour Lemrabet

An Impedance E¤ect of a Thin Adhesive Layer in Some Boundary Value and Transmission Problems Governed by Elliptic Di¤erential Equations

Keywords: Boundary value problem of elliptic type, Transmission problems, Impedance e¤ect, Thin layer. 2000 Mathematics Subject Classi…cation: 34K30, 35C15, 35J05, 35J25, 35Q60, 35Q61, 47D06. 1

In this paper we consider a problem of two bodies bonded through a thin adhesive layer (a third material) of thickness . Leting go to zero, one obtains a boundary value transmission problem set on a …xed domain. We then give new results for the study of this problem in the framework of Hölder spaces: an explicit representation of the solution and necessary and su¢ cient conditions at the interface for its optimal regularity are obtained using the semigroups theory and the interpolation spaces.

Introduction

Consider the boundary value and transmission second-order operational problem

(P ) 8 > > > > > > > < > > > > > > > : u 00 (x) + Au (x) = g (x) on ] 1; 0[ [ ]0; [ [ ] ; 1 + [ u ( 1) = f u 0 (1 + ) = f + u (0 ) = u (0 + ) + ; u = u + + p u 0 (0 ) = p 0 u 0 (0 + ) + a p 0 u 0 = p + u 0 + + b ; (1)
set in some complex Banach space E; here A is a closed linear operator of domain D(A) E (not necessarily dense in E) which veri…es the Krein's ellipticity (see Section 3,[START_REF] Tanabe | Equations of Evolution, Monographs and Studies in Mathematics 6[END_REF], f , f + , , , a , b are given in E and satisfy some necessary and su¢ cient conditions which will be speci…ed later.

The second member g is such that 8 > < > :

g = g j [ 1;0] 2 C ([ 1; 0] ; E) g 0 = g j [0; ] 2 C ([0; ] ; E) g + = g j [ ;1+ ] 2 C ([ ; 1 + ] ; E) ;
(with 0 < < 1). It is not di¢ cult to prove that the holderianity of g , g 0 and g + imply the global holderianity of g on [ 1; 1 + ] if and only if g (0) = g 0 (0) and g 0 ( ) = g + ( ):

We do not assume these two conditions. Set u = u j] 1;0[ , u 0 = u j]0; [ , u + = u j] ;1+ [ , then problem (P ) writes

(P ) 8 > > > > > > > > > > > > > > < > > > > > > > > > > > > > > :
(EQS) 8 > < > :

u 00 (x) + Au (x) = g (x) on ] 1; 0[ u 0 00 (x) + Au 0 (x) = g 0 (x) on ]0; [ u + 00 (x) + Au + (x) = g + (x) on ] ; 1 + [ ;

(BC) u ( 1) = f u + 0 (1 + ) = f + ; (T C) 8 > > < > > :
u (0) = u 0 (0) + u 0 ( ) = u + ( ) + p u 0 (0) = p 0 u 0 0 (0) + a p 0 u 0 0 ( ) = p + u + 0 ( ) + b :

The numerical solution of this problem is usually very di¢ cult to compute. In fact, the small thickness of the thin layer generates di¢ culties in the meshing. As ! 0, the interval ]0; [ degenerates into the point f0g and we can no more have an equation on it. The interval ] ; + 1[ becomes ]0; 1[ : Therefore, the main question is: what will be the appropriate transmission conditions at the interface point f0g which describe correctly the e¤ect of the thin adhesive layer ]0; [ as ! 0?

We will answer formally to this question in the most interesting case, characterizing so the e¤ect of the small bond ]0; [ (as ! 0).

Let us begin by giving a formal derivation of the e¤ect of the small bond ]0; [ (as ! 0).

In order to deal with our problem (P ), we solve the scalar equation on the small intervall ]0; [ and write down relations between the Cauchy data u 0 (0); u 0 0 (0) and u 0 ( ); u 0 0 ( ) :

Then making use of the transmission condition at f0g and f g, we will deduce relations linking u (0); p u 0 (0) ; to u + ( ); p + u + 0 ( ) ; which allow us, as ! 0, to obtain the limiting transmission conditions at f0g. Therefore, we will see that the interesting limit problem writes in the form

(P 1A ) 8 > > > > < > > > > : (u ) 00 (x) + Au (x) = g (x) on ] 1; 0[ (w + ) 00 (x) + Aw + (x) = h + (x) on ]0; 1[ u ( 1) = f ; (w + ) 0 (1) = f + u (0) = w + (0) + ' p (u ) 0 (0) p + (w + ) 0 (0) = qAu (0) + ;
see the details in Subsection 2.1.

Our main results concerning problem (P 1A ) are summarized in the following Theorem.

Theorem 1 Let g 2 C ([ 1; 0] ; E), h + 2 C ([0; 1] ; E) with 0 < < 1 and f 2 D (A), f + 2 D (( A)
1=2 ), ' 2 D(A), 2 E. Assume [START_REF] Tanabe | Equations of Evolution, Monographs and Studies in Mathematics 6[END_REF]. Then problem (P 1A ) has a unique solution

u (x) :] 1; 0[! E w + (x) :]0; 1[! E such that 1. u 2 C([ 1; 0]; D(A))\C 2 ([ 1; 0]; E), w + 2 C([0; 1]; D(A))\C 2 ([0; 1]; E) if and only if 8 > > > < > > > : ( s) Af g ( 1) 2 D(A) ( A) 1=2 f + 2 D(A) ( s) qg (0) + 2 D(A) h + (0) g (0) + A' 2 D(A): 2. Au (:); u 00 2 C ([ 1; 0] ; E) ; Aw + (:); w 00 + 2 C ([0; 1] ; E) if and only if 8 > > < > > : ( r) Af g ( 1) 2 D A ( =2; +1) ( A) 1=2 f + 2 D A ( =2; +1) ( r) qg (0) + 2 D A ( =2; +1) h + (0) g (0) + A' 2 D A ( =2; +1) :
In this main result, note that ( s) and ( s) are respectively the necessary and su¢ cient compatibility conditions at the boundary and the necessary and su¢ cient compatibility conditions at the interface f0g to obtain a strict solution u = (u ; w + ). Similarly, ( r) and ( r) are those to obtain optimal regularities on u:

The de…nition and the properties of the interpolation space D A ( =2; +1) are given, for instance, in ( [START_REF] Grisvard | Spazi di Tracce e Applicazioni[END_REF])

Many authors have worked on analogous problems, see [START_REF] Krasucki | Analysis of Interfaces of Variable Sti¤ ness[END_REF], [START_REF] Krasucki | Yield Design of Bonded Joints[END_REF] and [START_REF] Geymonat | Mathematical Analysis of a Bonded Joint with a Soft Thin Adhesive[END_REF] in hilbertian spaces. In [START_REF] Belhamiti | Transmission Problems in a Thin Layer Set in the Framework of the Hölder Spaces, Resolution and Impedance[END_REF] and [START_REF] Dore | An Abstract Transmission Problem in a Thin Layer, I: Sharp Estimates[END_REF], a study is given for a similar problem respectively in the framework of Hölder spaces and L p -spaces. These two last studies have considered only two materials. In our work we will use some techniques of these approachs which are based on the theory of semigroups, the Dunford functional calculus and the interpolation spaces.

This paper is organized as follows.

In Section 2, one gives the formal calculus for the limiting transmission problem and a concrete problem which motivates our study. In section 3, we give the basic hypothesis and some technical lemmas useful to the study of our problem (P 1A ). Section 4 is devoted to the derivation of an explicit representation of the solution of (P 1A ). In Sections 5 and 6 we study the solution and give in addition necessary and su¢ cient compatibility conditions on the data in order to obtain the above Theorem. In a last Section 7, we go back to the main physical example given in Subsection 2.2 and apply our results in the case of the space E = C 0 ( ) of the -Hölder continuous functions vanishing on the boundary @ : 2 Formal derivation of the limiting e¤ect of the thin junction

Derivation of the transmission conditions

In order to have an idea at least formally of the limiting problem, let us …rst consider the case when operator A is replaced by a complex scalar z (with z 2 CnR + ) and for the simplicity

g j]0; [ = 0 = = a = b = 0:
De…ne functions w + and h + on the …xed intervall ]0; 1[ by

w + (x) = u + ( + x); h + (x) = g + ( + x):
For simplicity, we have supposed that these functions do not depend on . The equation on the intervall ]0; [ writes u 0 00 + zu 0 = 0;

which gives u 0 (x) = C 1 e p zx + C 2 e p z( x) ;
where C 1 and C 2 are constants to be …xed by the boundary conditions. We thus have 8 > > > < > > > :

u 0 (0) = C 1 + C 2 e p z u 0 ( ) = C 1 e p z + C 2 u 0 0 (0) = C 1 p z + C 2 p ze p z u 0 0 ( ) = C 1 p ze p z + p zC 2 :
Now, the transmission conditions 8 > > < > > : The limiting transmission conditions as ! 0, which are obtained from the analysis of (2) and (3) depend on the behavior of p , p 0 and p + with respect to . So we must assume some conditions on p , p 0 and p + . The most interesting case we will consider is the following p and p + are independent of and p 0 = q ;

u (0) = u 0 (0) u 0 ( ) = u + ( ) = w + (0) p u 0 (0) = p 0 u 0 0 (0) p 0 u 0 0 ( ) = p + u + 0 ( ) = p + (w + ) 0 (0) ; lead to 2e p z p zw + (0) (2)
where q is a …xed positive number. This problem may model an electrostatic potential u in an heterogeneous material (see next subsection). The heterogeneity of the material is translated by the discontinuity of the conductivity coe¢ cient p. In this situation the small bond is highly conductive, implying thus the continuity of the potential through the sheet, but the normal component of the electric induction …eld is no longer continuous through the interface f0g, it has a jump proportional to the potential at f0g. In that case one has

1 e 2 p z 1 p 0 = 1 e 2 p z q ! 0 as ! 0, p 0 1 e 2 p z = q 1 e 2 p z ! 2q p z as ! 0,
and then the limiting transmission conditions are

w + (0) = u (0) p + (w + ) 0 (0) = q p z p zu (0) + p (u ) 0 (0);
and the limiting scalar problem becomes

(P 1z ) 8 > > > > < > > > > : (u ) 00 (x) + zu (x) = g (x) on ] 1; 0[ (w + ) 00 (x) + zw + (x) = h + (x) on ]0; 1[ u ( 1) = f ; (w + ) 0 (1) = f + u (0) = w + (0) p (u ) 0 (0) p + (w + ) 0 (0) = qzu (0): (5)
Remark 2 Another interesting case is when p and p + are constant with respect to but p 0 = q. In that case the small bond is not su¢ ciently conductive, implying thus a jump of the potential at f0g, but the normal component of the electric induction …eld remains continuous through the interface f0g and it is proportional to the jump of the potential at f0g. In that case, one has

1 e 2 p z 1 p 0 = 1 e 2 p z 1 q ! 2 1 q p z a s ! 0 , p 0 1 e 2 p z = q 1 e 2 p z ! 0 as ! 0;
then the transmission conditions are

qw + (0) = qu (0) + p + (u ) 0 (0) p + (w + ) 0 (0) = p (u ) 0 (0);
and the limiting problem becomes

(P 2A ) 8 > > < > > : (u ) 00 (x) + Au (x) = g (x) on ] 1; 0[ (w + ) 00 (x) + Aw + (x) = h + (x) on ]0; 1[ u ( 1) = f ; (w + ) 0 (1) = f + ; p (u ) 0 (0) = q (u (0) w + (0)) = p + (w + ) 0 (0)):
In the general case where the data g j]0; [ , , , a and b are di¤erent from 0, one obtains the following relations between w + (0); (w + ) 0 (0) and When we assume (4), we obtain the following limiting scalar problem with non homogeneous transmission conditions

u (0); u 0 (0) 2e p z p zw + (0) (6)
(P 1z ) 8 > > > > < > > > > : (u ) 00 (x) + zu (x) = g (x) on ] 1; 0[ (w + ) 00 (x) + zw + (x) = h + (x) on ]0; 1[ u ( 1) = f ; (w + ) 0 (1) = f + u (0) = w + (0) + ' p (u ) 0 (0) p + (w + ) 0 (0) = qzu (0) + : (8)
where Therefore, in this work we will focus ourselves on the complete analysis of the following problem

' = lim
(P 1A ) 8 > > > > < > > > > : (u ) 00 (x) + Au (x) = g (x) on ] 1; 0[ (w + ) 00 (x) + Aw + (x) = h + (x) on ]0; 1[ u ( 1) = f ; (w + ) 0 (1) = f + u (0) = w + (0) + ' p (u ) 0 (0) p + (w + ) 0 (0) = qAu (0) + : (9)

Electrostatic potential in a heterogeneous cylinder

r: pru = pg in G u = 0 on ] 1; 1 + [ u = f on f 1g p @u @x = f + on f1 + g ; (10)
models an electrostatic problem in G . The function u is the electrostatic potential, ru is the electric …eld and pru is the electric induction …eld. The heterogeneity of the material is translated by the discontinuity of the conductivity coe¢ cient p :

p = 8 < : p in ] 1; 0[ p 0 in ]0; [ p + in ] ; 1 + [ :
where p , p 0 and p + are positive constants. The function g is a given electric density, f is a …xed surface potential and f + a …xed surface induction. Set 8 > < > :

u = u j] 1;0[ ; g = g j] 1;0[ ; u 0 = u j]0; [ ; g 0 = g j]0; [ ; u + = u j] ;1+ [ , g + = g j] ;1+ [ ; then, the equation r: pru = pg in G ;
is equivalent to the following equations 8 > < > :

u = g in ] 1; 0[ u 0 = g 0 in ]0; [ u + = g + in ] ; 1 + [ ; (11)
with the transmission conditions

8 > > > > > > > < > > > > > > > : u = u 0 on f0g u 0 = u + on f g p @u @x = p 0 @u 0 @x on f0g p 0 @u 0 @x = p + @u + @x on f g : (12)
Let us de…ne the functions w + , h + on the …xed intervall ]0; 1[ by setting

w + (x; y) = u + ( + x; y); h + (x; y) = g + ( + x; y);
where we have assumed, for simplicity, that they do not depend on .

The approximations for u 0 (x; :) and @u 0 @x (x; :) as ! 0; give

u 0 ( ; y) ' u 0 (0; y) + @u 0 @x (0; y) + 2 2 @ 2 u 0 @x 2 (0; y) = u 0 (0; y) + @u 0 @x (0; y) 2 2
y u 0 + g 0 (0; y); @u 0 @x ( ; y) ' @u 0 @x (0; y) + @ 2 u 0 @x 2 (0; y) = @u 0 @x (0; y) y u 0 g 0 (0; y):

Using these relations with ( 11) and ( 12), we get the following problem set on the …xed domain (not depending on ):

(EQS) u = g in ] 1; 0[ w + = h + in ]0; 1[ ; (13)
under the boundary conditions

(BC) 8 > > > < > > > : u = 0 on ] 1; 0[ w + = 0 on ]0; 1[ u = f on f 1g p + @w + @x = f + on f1g ; (14)
and the following transmission conditions (depending on )

(T C) 8 > > <
> > :

w + = u + p p 0 @u @x 2 2 ( y u 0 g 0 ) on f0g p + @w + @x = p @u @x p 0 ( y u g 0 ) on f0g ;
which model the e¤ect of the thin cylinder ]0; [ on the other parts of the cylinder. The thin cylinder ]0; [ is modelled by the sheet f0g and its e¤ect is modelled through these transmission conditions.

There are two limiting cases of particular interest. The …rst case is

p 0 = q ;
(q is a positive constant) which is considered in this work assumes that the thin layer is highly conductive; as goes to 0, we obtain 8 <

:

w + = u on f0g p + @w + @x = p @u @x q ( y u + g 0 ) on f0g ; ( 15
)
which corresponds to the fact that the potential is continuous through the sheet, but the normal component of the electric induction …eld has a jump proportional to the potential.

The second case is p 0 = q and corresponds to the fact that the thin layer is poorly conductive. We get, as goes to 0: 8 > < > :

w + = u + 1 q p @u @x on f0g p + @w + @x = p @u @x on f0g ;
here, the normal component of the electric induction is continuous through the sheet but the potential has a jump proportional to the normal componant of the electric induction …eld. Therefore, using the classical operational notations u (x)(y) := u (x; y); w + (x) (y) := w + (x; y); :::

the concrete problem ( 13)-( 14)- [START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF] writes exactly in the form (9) with = qg 0 and ' = 0 and

D (A) = v 2 W 2;p (0; 1) : v (0) = v (1) = 0 Av (y) = v 00 (y) :
in the case E = L p (0; 1), or

D (A) = v 2 C 2 [0; 1] : v (0) = v (1) = 0 Av (y) = v 00 (y) : in the case E = C [0; 1] : 3

Hypotheses and technical lemmas

We assume in all this work the following ellipticity hypothesis: 8 < :

for any 2 ]0; [ (A) S [ f0g and
9C > 0 : 8 2 S [ f0g ( I A) 1 L(X) 6 C 1 + j j ; ( 16
)
where % (A) denotes the resolvent set of A and S = fz 2 Cn f0g : jarg zj < g :

This assumption implies that there exist a ball B (0; r 0 ), r 0 > 0, such that

% (A) B (0; r 0 ); (18)
and the estimate in ( 16) is still true in S [ B (0; r 0 ). It is well known that the above assumption implies that the square root p A is well de…ned and B = p A generates an analytic semigroup

e B >0 ; (19)
which is not necessarily continuous at 0. Let us recall the following important properties of B proved in [START_REF] Sinestrari | On the Abstract Cauchy Problem of Parabolic Type in Space of Continuous Functions[END_REF].

Lemma 3 Let 2 E and 2 C ([0; T ] ; E) with T > 0. Then 1. e B ! as ! 0 + i¤ 2 D(B) = D(A), (see [14], Proposition 1.2 ; p.20), 2. 7 ! e B 2 C ([0; T ] ; E) i¤ 2 D B ( ; +1) = D A ( =2; +1
), (see [START_REF] Sinestrari | On the Abstract Cauchy Problem of Parabolic Type in Space of Continuous Functions[END_REF], Proposition 1.12 ; p.29),

3. 7 ! R 0 Be ( t)B [ (t) ( )] dt 2 C ([0; T ] ; E) \ B (0; T ; D B ( ; +1)),
(see [START_REF] Sinestrari | On the Abstract Cauchy Problem of Parabolic Type in Space of Continuous Functions[END_REF], Theorem 4.5 ; p.53).

Let 2]0; =2[ and set S(r 0 ; ) = fz 2 Cn f0g : jzj > r 0 sin and jarg zj < g :

Note that for all w 2 S(r 0 ; ), one clearly has Re w > r 0 sin :

We will also use the following result proved in Proposition 4.10, p. 1880 in [START_REF] Dore | An Abstract Transmission Problem in a Thin Layer, I: Sharp Estimates[END_REF].

Lemma 4 For any w 2 S(r 0 ; ), one has 1. jarg (1 e w ) arg (1 + e w )j < ;

2. j1 + e w j > C = 1 e 2 tan > 0;

3. j1 e w j Re w 1 + Re w > r 0 sin 1 + r 0 sin = c 1 :

Now, consider the following space H 1 (S(r 0 ; )) = ff : f is an holomorphic and bounded function on S(r 0 ; )g , then, under our assumptions on A, if f 2 H 1 (S(r 0 ; )) is such that 1=f 2 H 1 (S(r 0 ; )) and (1=f )( B) 2 L(X), then f ( B) is invertible with a bounded inverse and

[f ( B)] 1 = (1=f )( B);
see, for instance [START_REF] Dore | H 1 Functional Calculus for Sectorial and Bisectorial Operators[END_REF] or [START_REF] Haase | The Functional Calculus for Sectorial Operators and Similarity Methods[END_REF], p. 45, Remark 2.5.1.

On the other hand we recall that operator I e 2B and I +e 2B are boundedly invertible; see for instance [START_REF] Lunardi | Analytic Semigroups and Optimal Regularity in Parabolic Problems[END_REF], p. 60, Proposition 2.3.6.

Let us now apply this result to the following operator for w 2 S(r 0 ; ) with Re w x .

= I p q B 1 I e 2B 1 I + e 2B p + q B 1 I + e 2B 1 I
In the compact sector K x = fw 2 Cn f0g : r 0 sin 6 Re w 6 x and jarg(w)j 6 g ; there is at most a …nite number of roots of f (w) (not belonging to [0; x ]), (see the remark on Proposition 4.1 , p. 41 in [START_REF] Cartan | Théorie Elémentaire des Fonctions Analytiques d'une ou Plusieurs Variables Complexes[END_REF]); so there exists 2]0; ] such that f (w) does not vanish on r0; := fw 2 Cn f0g : r 0 sin 6 Re w 6 x and jarg(w)j 6 g :

From ( 21), we conclude that f (w) does not vanish on S(r 0 ; ). Hence f 2 H 1 (S(r 0 ; )), 1=f 2 H 1 (S(r 0 ; )) and consequently (1=f )( B) 2 L(X) and

= f ( B) = I p q B 1 I e 2B 1 I + e 2B p + q B 1 I + e 2B 1 I e 2B
is boundedly invertible.

4 Representation of the solution of (P 1A )

It is well known that, the solution of the second order following equation

u 00 (x) B 2 u(x) = g(x); x 2]a; b[ in E; writes u(x) = e (x a)B + e (b x)B + v (g) (x) with ; 2 E and v (g) (x) = 1 2 x Z a e (x t)B B 1 g (t) dt + 1 2 b Z x e (t x)B B 1 g (t) dt: Therefore u (x) = e xB + e (1+x)B + v (g ) (x) ; x 2] 1; 0[ w + (x) = e xB + + e (1 x)B + + v (h + ) (x) ; x 2]0; 1[ with ; ; + ; + 2 E and v (g ) (x) = 1 2 x Z 1 e (x t)B B 1 g (t) dt + 1 2 0 Z x e (t x)B B 1 g (t) dt; v (h + ) (x) = 1 2 x Z 0 e (x t)B B 1 h + (t) dt + 1 2 1 Z x e (t x)B B 1 h + (t) dt:
The boundary conditions

u ( 1) = e B + + v (g ) ( 1) = f w 0 + (1) = Be B + B + + v 0 (h + ) (1) = f + ; give = e B v (g ) ( 1) + f + = e B + + B 1 v 0 (h + ) (1) B 1 f + ; (23)
and the transmission conditions

u (0) = w + (0) + ' p (u ) 0 (0) p + (w + ) 0 (0) = qAu (0) + ; give 8 > > > < > > > : + e B + e B + = ' v (g ) (0) + v (h + ) (0) p + e B p + + e B + + qB + e B = p B 1 v 0 (g ) (0) + p + B 1 v 0 (h + ) (0) qBv (g ) (0) + B 1 :
Using (23), one obtains the system

I e 2B

I + e 2B + = F p I + e 2B + qB I e 2B p + I e 2B + = G ;

(24)

where 8 > > > > > < > > > > > : F = ' e B f e B B 1 f + + e B v (g ) ( 1) v (g ) (0) + e B B 1 v 0 (h + ) (1) + v (h + ) (0) G = B 1 (p + qB) e B f + p + e B B 1 f + + (p + qB) e B v (g ) ( 1) p + e B B 1 v 0 (h + ) (1) p B 1 v 0 (g ) (0) + p + B 1 v 0 (h + ) (0) qBv (g ) (0) : (25)
The abstract determinant of this system (acting on D(B)) is

(B) (
1 p I + e 2B + qB I e 2B F I e 2B G +B 1 v 0 (h + ) (1) B 1 f + :
Therefore, the functions

u (x) = e xB + e (1+x)B + v (g ) (x) ; x 2] 1; 0[ w + (x) = e xB + + e (1 x)B + + v (h + ) (x) ; x 2]0; 1[; with v (g ) (x) = 1 2 x Z 1 e (x t)B B 1 g (t) dt + 1 2 0 Z x e (t x)B B 1 g (t) dt; and v (h + ) (x) = 1 2 x Z 0 e (x t)B B 1 h + (t) dt + 1 2 1 Z x e (t x)B B 1 h + (t) dt;
are completely determined.

5 Analysis of v (g ) et v (h + ) and their dérivatives

5.0.1 Analysis of v (g ) on ( 1; 0) One has v (g ) (x) = 1 2 x Z 1 e (x t)B B 1 [g (t) g (x)] dt + 1 2 0 Z x e (t x)B B 1 [g (t) g (x)] dt + 1 2 x Z 1 e (x t)B B 1 g (x) dt + 1 2 0 Z x e (t x)B B 1 g (x) dt = 1 2 x Z 1 e (x t)B B 1 [g (t) g (x)] dt + 1 2 0 Z x e (t x)B B 1 [g (t) g (x)] dt + 1 2 e (x+1)B B 2 g (x) + 1 2 e xB B 2 g (x) B 2 g (x) ; and v (g ) (x) = 1 2 x Z 1 e (x t)B B 1 [g (t) g (x)] dt + 1 2 0 Z x e (t x)B B 1 [g (t) g (x)] dt + 1 2 e (x+1)B B 2 [g (x) g ( 1)] + 1 2 e xB B 2 [g (x) g (0)] B 2 g (x) + 1 2 e (x+1)B B 2 g ( 1) + 1 2 e xB B 2 g (0) = v R (g ) (x) + v S; 1 (g ) (x) + v S;0 (g ) (x) ; where v S; 1 (g ) (x) = 1 2 e (x+1)B B 2 g ( 1) ; v S;0 (g ) (x) = 1 2 e xB B 2 g (0) : Due to Lemma 3, v R (g ) has the following maximal regularity properties ( v R (g ) 2 C([ 1; 0]; D B 2 ) B 2 v R (g ) 2 C ([ 1; 0]; E);
while for v S; 1 (g ) and v S;0 (g ) (x) we only have

( v S; 1 (g ) 2 C(] 1; 0] ; D B 2 ); B 2 v S; 1 (g ) 2 C (] 1; 0]; E) v S;0 (g ) 2 C([ 1; 0[; D B 2 ); B 2 v S;0 (g ) 2 C ([ 1; 0[; E):
The behaviour of B 2 v S; 1 (g ) (x) in the neighbourhood of 1 is that of

1 2 e (x+1)B g ( 1) ; (31)
and the behaviour of B 2 v S;0 (g ) (x) in the neighbourhood of 0 is that of

1 2 e (x+1)B g (0) : (32)
5.0.2 Analysis of v (h + ) on (0; 1)

Recall that v (h + ) (x) = 1 2 x Z 0 e (x t)B B 1 h + (t) dt + 1 2 1 Z x e (t x)B B 1 h + (t) dt;
which can be written as

v (h + ) (x) = 1 2 x Z 0 e (x t)B B 1 [h + (t) h + (x)] dt + 1 2 1 Z x e (t x)B B 1 [h + (t) h + (x)] dt + 1 2 e xB B 2 h + (x) + 1 2 e (1 x)B B 2 h + (x) B 2 h + (x) = 1 2 x Z 0 e (x t)B B 1 [h + (t) h + (x)] dt + 1 2 1 Z x e (t x)B B 1 [h + (t) h + (x)] dt B 2 h + (x) + 1 2 e xB B 2 [h + (x) h + (0)] + 1 2 e (1 x)B B 2 [h + (x) h + (1)] + 1 2 e xB B 2 h + (0) + 1 2 e (1 x)B B 2 h + (1) = v R (h + ) (x) + v S;0 (h + ) (x) + v S;1 (h + ) (x) with v S;0 (h + ) (x) = 1 2 e xB B 2 h + (0); v S;1 (h + ) (x) = 1 2 e (1 x)B B 2 h + (1): (33)
Due to Lemma (3) v R;0 (h + ) has the following maximal regularity properties

( v R (h + ) 2 C([0; 1]; D B 2 ) B 2 v R (h + ) 2 C ([0; 1]; E);
while for v S;0 (h + ) and v S;1 (h + ) we only have

( v S;0 (h + ) 2 C(]0; 1]; D B 2 ); B 2 v S;0 (h + ) 2 C (]0; 1]; E) v S;1 (h + ) 2 C([0; 1[; D B 2 ); B 2 v S;1 (h + ) 2 C ([0; 1[; E): Remark 6
We have in view the study of the regularity of u on [ 1; 0] and

w + on [0; 1]. Since u (x) = e xB + e (1+x)B + v (g ) (x) ; x 2] 1; 0[ w + (x) = e xB + + e (1 x)B + + v (h + ) (x) ; x 2]0; 1[
the singular parts v S; 1 (g ) ; v S;0 (g ) ; v S;0 (h + ) and v S;1 (h + ) must be associated with the other singular parts of u and w + to …nd necessary and su¢ cient conditions on the data '; ; f ; f + ; g and h + in order to get optimal regularity for u on [ 1; 0] and w + on [0; 1].

Analysis of the derivatives

In order to study the regularity of the solution of our problem, we also need to analyze the behaviour of the derivative of v (g ) and v (h + ). We have, for all

x 2 ] 1; 0[ v 0 (g ) (x) = 1 2 x Z 1 e (x t)B g (t) dt 1 2 0 Z x e (t x)B g (t) dt + 1 2 B 1 g (x) 1 2 B 1 g (x) = 1 2 x Z 1 e (x t)B g (t) dt 1 2 0 Z x e (t x)B g (t) dt; then v 0 (g ) (0) = 1 2 0 Z 1 e tB g (t) dt = 1 2 0 Z 1 e tB [g (t) g (0)] dt + 1 2 0 Z 1 e tB g (0) dt = 1 2 B 1 g (0) + 1 2 0 Z 1 e tB [g (t) g (0)] dt = 1 2 B 1 g (0) + R (g ) (0) :
Similarly, for all x 2 ]0; 1[, one has

v 0 (h + ) (x) = 1 2 x Z 0 e (x t)B h + (t) dt 1 2 1 Z x e (t x)B h + (t) dt; then v 0 (h + ) (0) = 1 2 1 Z 0 e tB h + (t) dt = 1 2 B 1 h + (0) 1 2 1 Z 0 e tB [h + (t) h + (0)] dt = 1 2 B 1 h + (0) + S (h + ) (0) ; and v 0 (h + ) (1) = 1 2 B 1 h + (1) + 1 2 1 Z 0 e (1 t)B [h + (t) h + (1)] dt (34) = 1 2 B 1 h + (1) + S (h + ) (1) :
6 Study of the regularity of the solution of problem (P 1A )

6.1 Necessary conditions on the transmission data

Assume that u and w + are strict solutions, that is (u )

0 ( 1 ) = lim 0 ! 1 u ( 1 ) u ( 0 ) 1 0 2 D(A);
we also deduce that (u ) 00 ( 1 ) 2 D(A):

Similarly, for any 2 2 [0; 1], (w

+ ) 0 ( 2 ); (w + ) 00 ( 2 ) 2 D(A): It follows that (u ) 00 (0) = g (0) Au (0) 2 D(A); (w + ) 00 (0) = h + (0) Aw + (0) 2 D(A):
Using the …rst transmission condition

' = u (0) w + (0) 2 D(A); one has (u ) 00 (0) = g (0) Au (0) = g (0) A' Aw + (0)
= g (0) A' h + (0) + (w + ) 00 (0); from which we deduce the …rst following necessary compatibility condition h + (0) g (0) + A' = (w + ) 00 (0) (u ) 00 (0) 2 D(A):

The second transmission condition gives

p (u ) 0 (0) p + (w + ) 0 (0) = + qAu (0) = + q [Au (0) g (0)] + qg (0) 2 D(A)
which implies the second following necessary compatibility condition

p (u ) 0 (0) p + (w + ) 0 (0) q [Au (0) g (0)] (36) = + qg (0) 2 D(A):
The two conditions (35)-( 36) are equivalent to

+ qg (0) 2 D(A) + q [h + (0) + A'] 2 D(A);
due to the identity

h + (0) g (0) + A' = 1 q [ + qh + (0) + qA'] 1 q [ + qg (0)] 2 D(A):
In the case = 0 and ' = 0, the necessary conditions on the transmission data are g (0) 2 D(A); h + (0) 2 D(A): Then the behaviour of B 2 u (x) near 1 is the same as that of the function

e (1+x)B B 2 [ v (g ) ( 1) + f ] + B 2 v (g ) (x) : One has v (g ) ( 1) = 1 2 0 Z 1 e (t+1)B B 1 g ( 1) dt + 1 2 0 Z 1 e (t+1)B B 1 [g (t) g ( 1)] dt = 1 2 B 2 g ( 1) + R(g ) ( 1)
;

and the term R(g ) ( 1) is regular since

B 2 [R(g ) ( 1)
[; E) i¤ g ( 1) + B 2 f 2 D(B 2 ) = D(B) B 2 u (:) 2 C ([ 1; 0[; E) i¤ g ( 1) + B 2 f 2 D B ( ; +1) :

Analysis of w + near 1

Let us assume the necessary condition f + 2 D (B). Recall that

w + (x) = e xB + + e (1 x)B + + v (h + ) (x) ; then the behaviour of B 2 w + (x) = B 2 h e xB + + e (1 x)B + + v (h + ) (x) i ;
near 1 is the same as that of

B 2 e (1 x)B + + v (h + ) (x) :
One has seen in ( 30), (34) that

+ = e B 1 p I + e 2B + qB I e 2B F + I e 2B G +B 1 v 0 (h + ) (1) B 1 f + ; with v 0 (h + ) (1) = 1 2 B 1 h + (1) + S (h + ) (1) : Since the behaviour of B 2 v (h + ) (x) near 1 is that of 1 2 e (1 x)B h + (1)
, as mentioned in (33), it follows that the behaviour of B 2 w + (x) near 1 is that of

e (1 x)B 1 2 h + (1) + Bf + + 1 2 e (1 x)B h + (1) = e (1 x)B Bf + :
Due to Lemma 3, B 2 w + (:) has the following maximal regularity near 1

( B 2 w + (:) 2 C(]0; 1]; E) i¤ f + 2 D (B) et Bf + 2 D (B) B 2 w + (:) 2 C (]0; 1]; E) i¤ Bf + 2 D B ( ; +1) :
6.4 Analysis of u at the interface 0

The behaviour of

B 2 u (x) = B 2 e xB + B 2 e (1+x)B + B 2 v (g ) (x) ;
near 0 is the same as that of

B 2 e xB + B 2 v (g ) (x) ;
and due to (32), it is the same as that

B 2 e xB + 1 2 e xB g (0) :
Now, in virtue of (27), one has as mentioned in the previous section, see (38). One then deduces that the behaviour of B 2 w + (x) near 0 is that of e xB qB 2 ' + qh + (0) :

= [ (B)] 1 p + I e 2B F I + e 2B G = [ (B)] 1 [p + F G ] e 2B [ (B)] 1 [p + F + G ] : Write (B) = p + I e
Then, using again [START_REF] Sinestrari | On the Abstract Cauchy Problem of Parabolic Type in Space of Continuous Functions[END_REF], the necessary and su¢ cient conditon to obtain the regularity of w + near 0 depends on the regularity of qB 2 ' + qh + (0) :

We conclude by the following Theorem . w + = u on f0g p @u @x p + @w + @x = q y u + qg 0 on f0g :

Here is a bounded domain of R n , n > 1; with a regular boundary :

In view to illustrate our abstract analysis, we are going, in this section, to explicit and interpret our impedance compatibility conditions Then A verifes ( 16) and D(A) coincides with the well known little Hölder continuous functions space h 0 ( ), see [START_REF] Campanato | Generation of Analytic Semigroups by Elliptic Operators of Second Order in Hölder Spaces[END_REF], p. 497. Let us point out that the boundary condition in space C 0 ( ) is essential. Otherwise the estimate in ( 16) is not veri…ed, see Example 3.1.33 given [START_REF] Lunardi | Analytic Semigroups and Optimal Regularity in Parabolic Problems[END_REF], p. 110.

One also has

  These two last relations link the Cauchy data w + (0); (w + ) 0 (0) (at the interface f0g) of the function w + (de…ned on [0; 1]) to the Cauchy data u (0); u 0 (0) (at the interface f0g) of the function u (de…ned on [ 1; 0]).

+

  g j]0; [ ; ; a ; b ; where ' g j]0; [ ; ;

  0; [ ; ; a ; b ; :

  Consider the cylinder G = ] 1; 1 + [ constituted by the junction of two homogeneous cylinders G = ] 1; 0[ and G + = ] ; 1 + [ bonded together by the thin cylinder G 0 = ]0; [ (here is a bounded domain of R n , n > 1; with a regular boundary ). Denote by (x; y) the generic variable in G . The transmission problem 8 > > > < > > > :

  26) = p + I e 2B I e 2B I + e 2B p I + e 2B + qB I e 2B = p + I e 2B 2 + p I + e 2B 2 qB I e 2B I + e 2B : One has (B) : D (B) ! E with (B) = p + I e 2B 2 p I + e 2B 2 + qB I e 2B I + e 2B = qB I e 2B I + e 2B I p q B 1 I e 2B 1 I + e 2B p + q B 1 I + e 2B 1 I e 2B : = qB I e 2B I + e 2B ; which is boundedly invertible by Lemma 4. It follows that = [ (B)] 1 p + I e 2B F I + e 2B G ; (27) + = [ (B)] 1 p I + e 2B + qB I e 2B F I e 2B G ; (28) and then = e B [ (B)] 1 p + I e 2B F I + e 2B G v (g ) ( 1) + f ; (29) + (30) = e B [ (B)]

u 2 C

 2 ([ 1; 0]; D(A))\C 2 ([ 1; 0]; E) and w + 2 C([0; 1]; D(A))\C 2 ([0; 1]; E); then ' = u (0) w + (0) 2 D(A); and for any 1 2 [ 1; 0]

6. 2

 2 Analysis of u (x) near 1 Assume the necessary condition f 2 D B 2 = D (A). Recall that u (x) = e xB + e (1+x)B + v (g ) (x) ; where, due to (23), one has = e B v (g ) ( 1) + f :

] 2 D

 2 B ( ; +1) ; see Lemma[START_REF] Cartan | Théorie Elémentaire des Fonctions Analytiques d'une ou Plusieurs Variables Complexes[END_REF]. Now from the study and the results on the behaviour of B 2 v (g ) (x) near 1, see (31), one concludes that the behaviour of B 2 u (x) near 1 is the same as that of e (1+x)B 1 2 g ( 1) + B 2 f + 1 2 e (x+1)B g ( 1) = e (1+x)B g ( 1) + B 2 f : Due to Lemma (3), B 2 u (:) has the following maximal regularity properties near 1 ( B 2 u (:) 2 C([ 1; 0

Theorem 8 2 .that 1 q qB 2 '

 822 Assume[START_REF] Tanabe | Equations of Evolution, Monographs and Studies in Mathematics 6[END_REF]. Let f + 2 D (( A) 1=2 ), ' 2 D(A) and 2 E and letg 2 C ([ 1; 0] ; E), h + 2 C ([0; 1] ; E) with 0 < < 1. Then 1. w + 2 C([0; 1]; D(A)) \ C 2 ([0; 1]; E) if and only if ( A) 1=2 f + 2 D(A) and qB 2 ' + qh + (0) 2 D(A); Aw + (:); (w + ) 00 2 C [0; 1] ; E) if and only if ( A) 1=2 f + 2 D A ( =2; +1) and qB 2 ' + qh + (0) 2 D A ( =2; +1) : By observing + qh + (0) 1 q [ + qg (0)] = h + (0) g (0) + A';one deduce the complete Theorem announced in the Introduction.7 Going back to the concrete exampleLet us go back to our concrete limiting problem(x;y) u = g in ] 1; 0[ (x;y) w + = h + in ]0; 1[ ; on ] 1; 0[ w + = 0 on ]0; 1[ u = f on f 1g u = f on f 1g

  qg (0) + 2 D(A) h + (0) g (0) + A' 2 D(A); qg (0) + 2 D A ( =2; +1) h + (0) g (0) + A' 2 D A ( =2; +1) ;in the case of the following vector valued Banach spaces C([ 1; 0]; C 0 ( )) and C([0; 1]; C 0 ( )); provided that g 2 C ([ 1; 0]; C 0 ( )) and h + 2 C([0; 1]; C 0 ( )) and all the compatibility boundary conditions are satis…ed. Consider the following operator de…ned in C 0 ( ) as ( D(A) = n v 2 C 2 ( ) : v; y v 2 C 0 ( ) o [Av] (y) = y v(y):

D(A) = D A ( 1 += v 2 Cv 2 C 2 n v 2

 12222 =2; +1) = v 2 C 2+ ( ) : v = y v = 0 on @ ; see[START_REF] Lunardi | Analytic Semigroups and Optimal Regularity in Parabolic Problems[END_REF], p. 110, Corollary 3.1.32 and Corollary 3.1.35. The interpolation D A (1+ =2; +1) is intended in the Banach space C 0 ( ). Therefore, in our subspace E = C 0 ( ) C 0 ( ); one hasD A ( =2; +1) = (D(A); C 0 ( )) 1 =2;+1 = (D(A); C 0 ( )) 1 =2;+1 2+ ( ) : v = v = 0 on @ ; C 0 ( ) 2+ ( ) : v = v = 0 on @ ], p. 31, Corollary 1.2.18. The conditions = qg 0 2 E, f 2 D (A) become ( qg 0 2 C 0 ( ) f C 2 ( ) : v; y v 2 C 0 ( ) o

  2B 2 + p I + e 2B 2 qB I e 2B I + e 2B (37) = p + I + p I qB + T ; D B 2 . Therefore, the behaviour of B 2 w + (x) near 0 is that ofB 2 e xB [ (B)]1 B 1 qB 2 ' + qh + (0) :

	where, clearly T is a very regularizing operator. Then, it su¢ ces to analyze
	B 2 e xB [ (B)]	1 (p + F		G ) +	1 2	B 2 (p + I + p I qB ) g (0) ;
	near 0: First, let us evaluate the term
	(p + F	G ) +	1 2	B 2 (p

+ I + p I qB ) g (0) ; where Q 2

by using the expressions in (24). One has

and

Then

Now, as ' 2 D B 2 and B 2 [(p + + p ) g (0) p + h + (0)] 2 D B 2 ; the behaviour of

The determinant

being invertible of inverse

(where 1 I + e 2B 1 I e 2B 1 is bounded and commutes with e xB ), one concludes that the behaviour of B 2 u (x) near 0 is that of

Then, due to [START_REF] Sinestrari | On the Abstract Cauchy Problem of Parabolic Type in Space of Continuous Functions[END_REF], the necessary and su¢ cient conditon to obtain the regularity of u near 0 depends of the regularity of the term

One gets the following Theorem.

We have used tha fact that

6.5 Study of w + at the interface 0

Recall that

where

(see ( 28)). Therefore, the behaviour of B 2 w + (x) near 0 is the same as that of

Now as (see (33)) the behaviour of B 2 v (h + ) (x) near 0 is that of 1 2 e xB h + (0) and (B) = p + I + p I qB + T ;

(see (??)), it su¢ ces to analyze the behaviour of

where as we have seen

+ S ; with regular terms R and S . Then

and Af g ( 1) ; qg (0) + ; qh + (0) + 2 D(A) signi…es 8 < :


Srdjan Stanković

Irena Orović
email: irenao@ac.me

Image Watermarking based on

Keywords: Digital image watermarking, space/spatial-frequency image analysis, Hermite functions, Hermite expansion coefficients

An image watermarking scheme that combines Hermite functions expansion and space/spatial-frequency analysis is proposed. In the first step, the Hermite functions expansion is employed to select busy regions for watermark embedding. In the second step, the space/spatialfrequency representation and Hermite functions expansion are combined to design the imperceptible watermark, using the host local frequency content. The Hermite expansion has been done by using the fast Hermite projection method. Recursive realization of Hermite functions significantly speeds up the algorithms for regions selection and watermark design. The watermark detection is performed within the space/spatial-frequency domain. The detection performance is increased due to the high information redundancy in that domain in comparison with the space or frequency domains, respectively. The performance of the proposed procedure has been tested experimentally for different watermark strengths i.e. for different values of the Peak Signal to

I. INTRODUCTION

The rapid development of digital technology brings a need to protect multimedia content from unauthorized usage. One of the common solutions is known as digital watermarking with applications ranging from digital rights management to integrity protection. It refers to specific information hiding based on embedding a secret signal into the multimedia data, under the constraints of imperceptibility, security and robustness to attacks [1]- [START_REF] Muharemagić | Survey of Watermarking Techniques and Applications[END_REF]. The owner of the protected content should be able to detect the watermark in order to prove ownership. An efficient digital watermarking procedure provides a good compromise between the watermark invisibility (which will preserve the host signal quality) and its robustness against common attacks (aiming to make the watermark undetectable).

Most of the existing image watermarking techniques can be categorized into spatial domain or transform domain algorithms [START_REF] Barni | Watermarking Systems Engineering[END_REF]- [START_REF] Kwitt | Efficient Detection of Additive DWT Watermarking in the DWT domain[END_REF]. Watermark embedding in transform domains, e.g. frequency domains (through the Discrete Cosine Transform (DCT), Discrete Fourier Transform (DFT), Discrete Wavelet Transform (DWT), etc.), generally increases the robustness to attacks [START_REF] Nikolaidis | Asymptotically optimal detection for additive watermarking in the DCT and DWT domains[END_REF]- [START_REF] Nafornita | A New Pixel-Wise Mask for Watermarking[END_REF] compared to the spatial domain. However, the frequency domain embedding usually causes an uncontrolled spread of the watermark over a possibly large spatial region. It means that the watermark embedded in some components in one domain will be spread over all components in other domain. In this sense, special masking algorithms are required to preserve the watermark imperceptibility. In order to exploit the advantage of using both domains simultaneously, the joint space/spatial-frequency based watermarking has been considered [START_REF] Mobaseri | Digital watermarking in the joint time-frequency domain[END_REF]- [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF]. The space/spatial-frequency analysis provides a local frequency content estimation and analysis for each individual image pixel. Hence, the watermark spatial location and its frequency characteristics are controlled simultaneously. Furthermore, the space/spatial-frequency based procedure provides a high robustness to attacks and good detection performance. For each pixel in the spatial domain, we produce a block of space/spatial-frequency coefficients whose size depends on the window used for the calculation of space/spatial-frequency representation. It means that the number of coefficients used in detection is significantly increased comparing to the number of watermarked pixels: for each pixel, we have approximately 200 space/spatialfrequency coefficients for correlation (50% of the 20x20 blocks used in our case). Consequently, a successful detection in space/spatial-frequency domain usually requires smaller number of watermarked pixels than the procedures in DCT and DWT domain.

This work represents a modified and improved version of the space/spatial-frequency based image watermarking procedure [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF]. Prior to watermark embedding, a set of suitable pixels is determined. Generally, the imperceptibility-robustness compromise is significantly improved by embedding the watermark within the busy image regions. According to the masking effects of the human visual system (HVS), the busy or high contrast image areas can mask any distortion better than smooth areas, and thus, can support higher watermark signal without visibility problem [START_REF] Podilchuk | Image adaptive watermarking using visual models[END_REF]- [START_REF] Brett | Digital Watermarking using Improved Human Visual System Model[END_REF]. The HVS frequency sensitivity and masking are closely related, such that both effects suggest that we are less sensitive to noise in high frequencies producing the busy parts of the image [START_REF] Brett | Digital Watermarking using Improved Human Visual System Model[END_REF], [START_REF] Bradley | A wavelet visible difference predictor[END_REF]. The similar design [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF] uses the energy of local space/spatial-frequency content to determine the occupancy and busyness of the region. This requires exhaustive procedure for setting the thresholds, since the energy within the region is not exactly a natural and precise measure to characterize region busyness. Therefore, we propose a simple and efficient method to detect busy (high-variance) regions that is based on an incomplete Hermite functions expansion [START_REF] Krylov | Fast Hermite projection method[END_REF]. Each image region is firstly reconstructed by using a few low-order Hermite functions, and the mean square error (MSE) between the original and reconstructed region is used as a distinguishability criterion. The low-order functions are related to the low-frequency content, while the higher order functions have more high-frequency components [START_REF] Yuan | Conditions for Generation of Stable and Accurate Hybrid TD-FD MoM Solutions[END_REF], [START_REF] Brown | Exoplanets or Dynamic Atmospheres? The Radial Velocity and Line Shape Variations of 51 Pegasi and tau Bootis[END_REF], [START_REF] Krylov | Fast Hermite projection method[END_REF], [START_REF] Kortchagine | Projection filtering in image processing[END_REF].

Hence, when using a few low-order functions for the reconstruction, the MSE for busy regions will be pretty high, since these regions are characterized by high-frequency content [START_REF] Won | Stochastic Image Processing (Information Technology: Transmission, Processing and Storage)[END_REF]. More details the regions have, more functions are required for proper reconstruction.

The selected busy regions are used for watermark design and embedding. In order to identify space/spatial-frequency components that are, from the perceptual point of view, appropriate for watermark embedding, the Hermite projection method is applied in the space/spatial-frequency domain, as well. The most significant components (from the perceptual point of view) are omitted to avoid image quality degradation, and the watermark is designed according to the remaining components. The space/spatial-frequency content of the watermark changes for each selected pixel according to its local space/spatial-frequency content. This further means that the space/spatial-frequency content of the watermark is not uniform across the image, as it is the case in [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF], where the watermark is fixed to the middle frequency ring. Comparing to [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF], the proposed pixels selection and watermark modeling techniques improves the quality of watermarked image, which has been measured by the local PSNR.

The watermark detection is performed in the space/spatial-frequency domain, since it provides a large number of elements available for the correlation based detection. In this case, the information about watermark will be spread over a large number of coefficients in the space/spatial-frequency domain, even when a small number of pixels is used. In other words, the detection is flexible regarding the number of considered pixels, because even for a small number of pixels there is still large number of elements in correlation. The proposed procedure is tested on various images in the presence of different commonly used attacks.

The paper is organized as follows. The procedure for busy regions selection is given in Section II. The watermarking procedure that includes watermark design, embedding and detection is proposed in Section III. The efficiency of the proposed procedure is tested experimentally in Section IV, while the concluding remarks are given in Section V.

II. SELECTION OF WATERMARKING REGIONS

A. Hermite Projection Method

The Hermite projection method has been introduced in various image and speech processing applications [START_REF] Krylov | Fast Hermite projection method[END_REF]- [START_REF] Kortchagine | Projection filtering in image processing[END_REF]. The Hermite functions are used as the basis functions due to their favorable properties (orthogonality and good localization property [START_REF] Krylov | Fast Hermite projection method[END_REF], [START_REF] Eberlein | A new method for numerical evaluation of the Fourier transform[END_REF]), providing a unique representation of different signals. They provide better computational localization in both signal and transform domains in comparison with the trigonometric functions. Also, as it will be shown later, the coefficients of expansion are easily computed, especially having in mind the recursive procedure for Hermite functions calculation. Furthermore, by using the Hermite expansion, the signal energy is approximated by the numerical integral of the Gauss-Hermite type and converges more rapidly than the rectangle rule involved in the DCT computation [START_REF] Lazaridis | Signal compression method for biomedical image using the discrete orthogonal Gauss-Hermite transform[END_REF].

The Hermite functions have been defined as: .

(

) 1
To speed up the calculation procedure, they can be defined recursively as follows [START_REF] Kortchagine | Projection filtering in image processing[END_REF]:

Ψ = Ψ = - Ψ = Ψ - Ψ ∀ ≥ (2)
For the purpose of image analysis, the expansion into Hermite functions can be applied to image rows, columns or both. For the sake of simplicity, we consider the Hermite expansion of image rows (one row at a time). Therefore, an image is represented by using a given number K of Hermite functions as follows [START_REF] Kortchagine | Projection filtering in image processing[END_REF]:

where I(x,y) denotes an image of size P×Q (or a given region), while c p are coefficients of Hermite expansion. Note that, prior to the Hermite expansion, removing the baseline from I(x,y)

is required, since ( ) 0, p x x ψ →
→∞ for all p. The baseline is defined as:

( , ) ( ,0) ( ) ( , 0) x I x Q I x B y I x y Q - = + ⋅ . (4)
A complete image reconstruction can be achieved when the number of Hermite functions is equal to P. However, it has been shown that the image can be successfully reconstructed even by using a smaller number (K<P) of Hermite functions. Additional savings can be achieved by applying the Hermite projection method over both coordinates. The Hermite expansion coefficients can be calculated by using the Gauss-Hermite quadrature, as follows:

2 1 1 ( ) 1 ( ), ( ( ))
where m x are zeros of the Hermite polynomials, while I y (x m ) is I(x m ,y) for a fixed y.

B. A procedure for busy regions selection based on the Hermite Projection Method

Image analysis is usually based on texture feature extraction methods, and the most significant among them are the auto-correlation, edge frequency, primitive-length, Law's method and cooccurrence matrices [START_REF] Sharma | Evaluation of Texture Methods for Image Analysis[END_REF], [START_REF] Clausi | An analysis of co-occurrence texture statistics as a function of grey level quantization[END_REF]. Each of these methods uses different features for texture classification, while their performance depends on the quality of feature and amount of training data. The standard texture classification methods are mainly based on matching algorithms with a predefined texture patterns, patches or primitives [START_REF] Sharma | Evaluation of Texture Methods for Image Analysis[END_REF], [START_REF] Kandaswamy | Efficient Texture Analysis of SAR Imagery[END_REF]. The classification is then performed by determining the well-known co-occurrence matrices features [START_REF] Haddon | Co-occurrence matrices for image analysis[END_REF] (directionality, coarseness, granularity, etc.), which are compared between the observed image regions and patterns. The first limitation that we face with in our application is that it is hard to define a pattern or patch for busy regions in natural images, due to their diversity. Furthermore, the statistical measures extracted from co-occurrence matrices need to be combined to provide regions classification, while the thresholds associated with each measure are not standardized and they are subject of exhaustive experimentation.

In watermarking application we are interested in determining busy regions with large variations and without fine structures, or busy regions where different texture elements intersect. Moreover, we should also analyze the amount of variations within the region to decide whether it is suitable for watermarking. If we consider histogram analysis to characterize the variations within the region, we lack the information about the spatial distribution of different values. When observing the histograms of different regions (as shown later in Section IV), we may say that the histogram does not provide such obvious classification features. Moreover, we would need to perform detailed histogram analysis and certainly to define a completely new procedure for region classification.

The proposed approach provides a simple and reliable procedure that identifies the regions with fast variations and irregular structures, by using image independent threshold. The MSE between the original region and its reconstructed version provides a natural, simple and efficient measure of region complexity. More complex and disordered region produces higher MSE.

Moreover, it has been shown that the gap between flat (flat and moderately flat) and busy (busy and moderately busy) regions is significant. It is important to note that the proposed classification method has not been used so far, although the Hermite projection method has been applied in texture classification, but with a full set of functions and combined with additional procedures.

An image can be decomposed into a set of regions/blocks R i (x,y) of size P×Q. Each region can be separately decomposed into Hermite functions. When the region is uniform/flat, it can be reconstructed by using much fewer functions than for the busy regions. Hence, when the same number of Hermite functions is used for all image regions, the mean square error between the original and reconstructed region will be smaller if the region is more uniform. This property can be used as a criterion for an efficient busy region selection. The mean square error between the original and reconstructed region is given by:

2 1 0 1 { ( , )} ( , ) ( ) K i i p p x y p MSE R x y R x y c x PQ ψ - =     = -     ∑∑ ∑ , ( 6
)
where K is the number of Hermite functions used for reconstruction. Generally, the region size is firstly chosen approximately to correspond to the size used in texture features characterization (e.g., between 15x15 and 30x30). Then, based on a number of experiments, the optimal performance of the proposed procedure is achieved by using the regions 20x20. The values of K are determined experimentally by using a set of test images. Namely, for the observed 20x20 regions we use a fixed number of five Hermite functions and coefficients, which is 25% of number of pixels per line. This value of K is enough to reconstruct flat regions with small error and insufficient to reconstruct busy regions. Apparently, when working with such regions the classification performance does not vary significantly for different images.

The considered regions are classified as follows:

{ ( , )} ( , ) { ( , )} ( , ) i i i i MSE R x y R x y busy region MSE R x y R x y flat region >Ξ ⇒ ∈ ≤Ξ ⇒ ∈ ( 7
)
where Ξ is some predetermined threshold value (which will be discussed in Section IV).

Generally, the pixels from busy regions are more suitable for watermarking, since the changes in these regions are more difficult to perceive than in the flat regions. Consequently, watermark embedded in the busy regions can be stronger, and it is more likely to be robust to different attacks. In the proposed approach, for watermark embedding we use only the pixels from the busy regions.

III. WATERMARKING PROCEDURE

A block diagram of the proposed watermarking procedure is shown in Fig 1 .a. The image is divided into regions processed by using the Hermite expansion method to select busy image regions. The selected regions are used as input to the space/spatial-frequency processing block together with the original image to perform watermark modeling and embedding. When choosing the size of image regions, two important issues are considered: 1. small image regions are difficult for characterization in both spatial and frequency domain; 2. large regions can not be associated with local frequency content of pixels anymore. Hence, as discussed in the previous section, the optimal size is set to 20x20.

The space/spatial-frequency processing block is shown in Fig 1 .b. It is important to emphasize that the block scheme given in Fig. 1, includes two major steps from the calculation standpoint:

1. Hermite expansion method, 2. space/spatial-frequency transform and its inverse. The fast Hermite projection method [START_REF] Krylov | Fast Hermite projection method[END_REF], suitable for hardware realization, has been used for Hermite expansion. Although, the Hermite functions can be calculated recursively, we can even store and load functions from registers since we use only a first few Hermite functions. Concerning the space/spatial-frequency transform and its inverse, it is important to emphasize that the real-time hardware implementation for space/spatial-frequency image analysis (including the spacevarying filtering) has been already done in [START_REF] Stankovic | System architecture for space-frequency image analysis[END_REF], [START_REF] Stankovic | An architecture for the VLSI design of systems for time-frequency analysis and time-varying filtering[END_REF]. Additional operations required by the proposed procedure are image decomposition into 20x20 blocks and additive watermark embedding, with almost negligible computational complexity.

A. Watermark design and embedding

After the pixels that belong to the busy image regions are selected, the watermark embedding procedure is performed in the space/spatial-frequency domain. As a suitable space/spatialfrequency representation, the two-dimensional (2D) short-time Fourier transform (STFT) is considered. By using the 2D STFT, a region of pixels, captured by the window around the central one, will be mapped into the 2D frequency representation. Hence, to provide a certain model for watermark characteristics, for each pixel we use a 2D representation of its local frequency content. The watermark embedding will be done in the space/spatial-frequency domain, as well.

The 2D form of the STFT for an image I(x,y) is defined as follows:

2 / ( , ) ( ) ( ) ,

j k N I k STFT n h k I n k e π ω ω - = + ∑ ( 8
)
where h is a sliding window (N×N). The following notation is used:

1 2 ( , ) n n n = , { } 1 2 1 2 ( , ) , [ / 2, / 2 1] k k k k k N N = ∈ - -, 1 1 2 2 ( ) ( , ) n k n k n k + = + + , 1 1 2 2 ( ) k k k ω ω ω = + .
The STFT is calculated for a region around each pixel that was previously selected for watermarking. Thus, the windowed part of signal, centered at the selected pixel, is transformed into the frequency domain. The local frequency content is further used to provide a suitable mask, which will shape the watermark, as it is described below. In that sense, we should avoid the high-energy space/spatialfrequency image content. For this purpose, the Hermite expansion is employed within the space/spatial-frequency domain. By using a small number K q of Hermite functions, we may characterize the high energy spectral content. Then the remaining space/spatial-frequency content is obtained as follows: 1 2 1 2

1 2 1 2 1 2 1 2 0 ( , , , ) ( , , , ) ( ), K q I q q q STFTx n n STFT n n c for fixed ω ω ω ω ψ ω ω - = = = -∑ . ( 9
)
After the high energy components are removed, the space/spatial-frequency content 2 ( , )

STFTx n ω is used to create a space/spatial-frequency mask:

2 1, ( , ) 0 ( ) 0, . w STFTx n L otherwise ω ω  >  =    (10)
Then, the space/spatial-frequency characteristics of a watermark are modeled as follows: An illustration of the space/spatial-frequency mask and the corresponding space/spatialfrequency characteristics of the watermark are shown in Fig. 2, for one pixel from the busy image region (region with majority of busy pixels).

( , ) ( , ) ( , )
m w w w STFT n L n STFT n ω ω ω = , (11) i.e. ( , ) ( , )
Additive watermark embedding procedure is performed in the space/spatial frequency domain, as follows:

( , ) ( , ) ( , )

m Iw I w STFT n STFT n STFT n ω ω α ω = + , ( 12
)
where Note that in order to retrieve the signal from its space/spatial-frequency representation, the concept of space-varying filtering is applied. This technique has been introduced for the purpose of space-varying filtering of 2D signals [START_REF] Stankovic | On the local Frequency, Group Shift and Cross Terms in Some Multidimensional Time-Frequency Distributions; A Method for Multidimensional Time-Frequency Analysis[END_REF], [START_REF] Lj | Space/Spatial-Frequency Based Filtering[END_REF]. The input of the space-varying filter is Iw STFT , while the support function is defined as:

2 1 2 1 2 1, ( , ) 0 ( , ) . 0, I if STFT L otherwise ω ω ω ω  > =   (13)
Finally, the watermarked image is obtained as follows:

( ) ( , ) ( , )

w Iw I n STFT n L n ω ω ω = ∑ . ( 14
)

B. Watermark detection

The watermark detection is performed, within the space/spatial-frequency domain, by using the space/spatial-frequency domain correlation as in [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF]:

{ } ( , ) ( , )
( , ) ( , ) ( , )

w w m I w n I w w n Det STFT n STFT n STFT n L n STFT n ω ω ω ω ω ω ω ∑ ∑ ∑ ∑   = =       =     . ( 15
)
Note that the two-dimensional space/spatial-frequency region is centered on each selected pixel. Instead of a single pixel in the spatial domain, we use a STFT matrix with a certain number of elements that depends on the window size. Then the actual number of components used in the detection is much larger, providing a significant detection improvement. However, note that ( , )

m w
STFT n ω depends on the support function ( , ) w L n ω , which is specific for each pixel. Thus, in order to use [START_REF] Briassouli | Locally optimum nonlinearities for DCT watermark detection[END_REF], we would need to store the 2D watermark for each watermarked pixel. In other words, the 2D support function should be available for each pixel as side information. Hence, to provide a blind detection and to avoid large side information, which is usually undesirable in watermarking, we use the modified detector form that does not need any knowledge about the 2D support functions. The introduced modification does not affect significantly the detection results.

Hence, the detector form includes only the STFT of the starting watermark, created as a 2D pseudo-random sequence, which is the same for all pixels:

( , ) ( , )

w I w n Det STFT n STFT n ω ω ω ∑ ∑   =    
, and does not require knowledge about the support function. The detection performance is tested by using the following measure of detection quality [START_REF] Heeger | Signal Detection Theory[END_REF], [START_REF] Wickens | Elementary Signal Detection Theory[END_REF]:

where Det and σ 2 represent the mean value and the variance of the detector responses, w and z indicate watermarks (right keys or trials) and wrong trials, respectively. The probability of error P err can be easily calculated by using the measure MD (given by ( 16)), as follows:

1 ( ). 2 2 err MD P erfc = (17)
The normal distribution of detector's responses is assumed. Recall that the watermark is created from a pseudo-random sequence and for such watermarks the correlation based detector form can be safely assumed to be sufficiently independent [START_REF] Tefas | Performance analysis of correlation-based watermarking schemes employing Markov chaotic sequences[END_REF]. According to the central limit theorem, if n random variables distributed according to any probability density function (pdf) are combined, then their sum will have pdf which tends to a Gaussian for sufficiently large n [START_REF] Tefas | Performance analysis of correlation-based watermarking schemes employing Markov chaotic sequences[END_REF], which is the case we have. Probability of detection is obtained as P d =1-P err . By increasing the value of measure MD the probability of detection error decreases.

IV. EXPERIMENTAL RESULTS

The performance of the proposed procedure is tested experimentally using various standard test images, all of size 256x256. According to the practice in the literature (e.g., [START_REF] Nikolaidis | Asymptotically optimal detection for additive watermarking in the DCT and DWT domains[END_REF], [START_REF] Briassouli | Locally optimum nonlinearities for DCT watermark detection[END_REF], [START_REF] Ali Akhaee | Contourlet-Based Image Watermarking Using Optimum Detector in a Noisy Environment[END_REF], [START_REF] Akhaee | Robust Scaling-Based Image Watermarking Using Maximum-Likelihood Decoder With Optimum Strength Factor[END_REF]) we have tested the proposed procedure to a certain categories of test images.

Since the watermarking procedure aims to protect natural images, we have selected some specific and commonly used examples (most of them are standard Matlab test images) such as portraits (e.g. Lena, Elaine, Barbara), Landscapes (e.g., Boat, Lake, Bridge, Autumn), images with plenty of details or repeating shapes (e.g., Baboon, Tiger, Peppers, Pears), images with dominant brightness (F16).

Numerous realizations with different watermarks (right keys) and wrong trials have been done.

The STFT is calculated by using the regions of size 20×20. The pixels belonging to busy image regions are selected by using the proposed procedure based on the Hermite projection method.

The regions classification is done by using parameter 40 Ξ= . Note that the classification performance is not so sensitive to the selection of Ξ . Namely, the experiments have shown that the MSE (given by ( 6)) for flat and busy regions have significantly different values, and thus we can set parameter Ξ easily and flexibly for all tested images. In Fig. 3 we illustrate four different types of regions classified according to the value of MSE:

-Strictly flat regions have very low MSE (mostly <10), -Moderately flat regions may contain certain peripheral changes and thus slightly higher MSE, -Moderately busy regions contain quite noticeable dynamics (MSE>100), -Busy regions contain significant dynamics (MSE>200).

Note that there is a large difference between the MSE of moderately flat and moderately busy regions, which facilitates setting of the threshold value Ξ .

Additionally, we may observe that it would be difficult to conclude anything about regions dynamics on the basis of their histograms, since they do not have any visible dominant characteristics for a certain region category (Fig. 3). The proposed method for pixel selection provides m candidates for watermark embedding, but the embedding procedure can be done in any set of n points, where n≤m. Hence, we may fix the number of watermarked pixels for both low and high resolution images. The selected pixels (black points) for some sample images are given in Fig 4 . The watermark has been modeled according to [START_REF] Sattar | A new time-frequency based private fragile watermarking scheme for image authentication[END_REF] and embedded in the selected pixels within busy image regions.

The watermarking procedure has been performed for a certain number of watermarks (right keys) and wrong trials (wrong keys). In Fig. 5, we present the values of measure MD obtained for 40 right keys and 400 wrong trials. Note that the blind detection is used, i.e. the original image is not required in the detection process. The original and watermarked images Lena and Elaine are shown in Fig. 6. Furthermore, the procedure is tested in the presence of various attacks with different strengths.

The watermark embedding is done with a PSNR ranging in 47-50 dB. As an additional measure of image quality, we have calculated the Structural Similarity Measure (SSIM) [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF], [41] between the original and watermarked images, and the achieved value for different test images is ≥0.99 (out of 1), which means almost negligible difference between the watermarked and original image.

The procedure is tested under JPEG compression with various quality factors (lower quality factor (QF) corresponds to higher compression), impulse noise with different values of density parameter (0.0025, 0.005, 0.01, 0.05), Gaussian noise with different values of variance (0.001, 0.005, 0.01, 0.02). For example, the robustness to Gaussian noise is achievable for variance up to 0.02 (the reference value for noise variance is usually 0.01 [START_REF] Rahmani | A New Robust Watermarking Scheme to Increase Image Security[END_REF]- [START_REF] Li | Reliable Adaptive Watermarking Scheme Integrated with JPEG[END_REF]). Further increasing of variance introduces serious image quality degradations and makes an image unreadable and unusable. Furthermore, the robustness is tested in the presence of mixed Gaussian and impulse noise, median filtering, combined noise and filtering attacks, image darkening and lightening, cropping of rows and columns, image resizing with interpolation to the nearest neighbour, and a few combined attacks.

The proposed approach provides low probabilities of error for most of the tested attacks. On the other hand, it is worth to say that the proposed procedure preserves high image quality and provides high PSNR (approximately 50 dB), which is larger than in the other robust image watermarking algorithms [START_REF] Kumsawat | A Robust Image Watermarking Scheme Using Multiwavelet Tree[END_REF], [START_REF] Nafornita | A New Pixel-Wise Mask for Watermarking[END_REF], [START_REF] Kang | Improving robustness of quantization-based image watermarking via adaptive receiver[END_REF], [START_REF] Ali Akhaee | Contourlet-Based Image Watermarking Using Optimum Detector in a Noisy Environment[END_REF]. Nevertheless, the image quality will be compared later with the similar approach in terms of the local quality measure. Note that the proposed approach provides robustness to cropping of image rows and columns, as well as robustness to certain scaling (17%) and interpolation, with still low probability of error. However, the procedure is not robust to image rotation. This would require combination with certain templates that would be invariant to rotation. Lightening 40% 10 -14 10 -14 10 -7 10 -5 10 -6

A. Comparison results

The proposed procedure is compared with the previous space/spatial-frequency (SSF) based approach given in [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF], with the standard DCT based watermark embedding procedure, and DWT approaches based on the Generalized Gaussian and Cauchy pdf model. The results are reported in Table 2. The watermark is embedded with the same PSNR≈50dB in all cases. Some details for these procedures are given in the sequel.

1) The method given in [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF] assumes different region selection and watermark embedding regions compared to the proposed approach. The results are given in the Table 2 (3 rd column).

Although, some of the results achieved by the proposed method are quite similar to the results in [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF], the proposed method is also robust to some geometric attacks unlike [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF]. Additional advantages over the method [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF] will be discussed in Section IV.B.

2) The DCT based image watermarking procedure is considered as follows. According to [START_REF] Briassouli | Locally optimum nonlinearities for DCT watermark detection[END_REF], [START_REF] Cheng | An additive approach to transform-domain information hiding and optimum detection structure[END_REF] all low-to-middle frequency DCT coefficients are used for watermarking (22050 coefficients for images of size 256×256, which is approximately four times the number of watermarked pixels in the case of the proposed procedure). The detection is performed by using the traditional correlation detector. The results are given in the Table 2, 4 th column.

3) Finally, we consider the standard DWT domain watermarking. The watermark detection in the DWT subbands is performed by using detectors based on the Generalized Gaussian or Cauchy pdf [START_REF] Nikolaidis | Asymptotically optimal detection for additive watermarking in the DCT and DWT domains[END_REF], [START_REF] Kwitt | Efficient Detection of Additive DWT Watermarking in the DWT domain[END_REF], since in this case they provide better performance compared to the standard correlator. Namely, it is commonly assumed that the watermarked coefficients' pdf may correspond either to Generalized Gaussian or Cauchy function, and thus detection is performed by using one of the mentioned models [START_REF] Nikolaidis | Asymptotically optimal detection for additive watermarking in the DCT and DWT domains[END_REF]. A number of 4000 DWT coefficients has been used.

Further increasing of the coefficients number will not influence significantly the detection measure. The results are given in Table 2 (5 th and 6 th column).

Note that the detection results achieved by the proposed procedure outperform the DCT and DWT based procedure. Furthermore, we provide a comparison with [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF] in terms of the perceptual quality.

B. Comparison with the previous space/spatial-frequency approach [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF]:

1) Improved imperceptibility: Comparing to the previous space/spatial-frequency based image watermarking, the proposed approach provides better watermark imperceptibility, measured by the local PSNR (LPSNR). The time-frequency mask obtained by using Hermite functions expansion allows us to create a watermark that fits well into the corresponding space/spatialfrequency content of the local image region. In other words the watermark is better suited to the host signal content than in the case when the constant middle frequency content defined by the circular ring:

{ } 1 2 1 1 2 2 ( , ) : , n m D ω ω ρ ω ω ρ = < <
is used, as in [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF].

In order to prove the above considerations, we have calculated the local PSNR for each watermarked region, in the case of the proposed and the previous approach [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF]. Note that this local metrics provides much reliable description of image quality than the global PSNR.

For a fair comparison, in these experiments we have used the same watermark for both approaches, providing the same detection performances. Also, the same image regions have been used for watermarking. Then, the LPSNR is calculated and compared for watermarked image regions. The results have shown that the proposed approach provides higher local PSNR for approximately 90% of watermarked pixels/regions. In Table 3, we have illustrated the space/spatial-frequency content for a few randomly selected pixels. The values of LPSNR are reported as well, for the proposed approach and the approach in [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF]. The pixels within 20×20 region are used to obtain local space/spatial-frequency content of the central pixel (the surrounding pixels have similar space/spatial-frequency representation). Similar results are obtained for other pixels/regions.

2) Robustness to Sensing: Unlike in [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF], where the watermark is fixed to the middle frequency ring, in the proposed approach the space/spatial-frequency content of the watermark changes over different regions and adapts to the local space/spatial-frequency content of host region. It means that watermark estimation as an attempt to remove it (sensing attack), is practically impossible without knowing the parameters K and K q in (6) and ( 9), which are used as security keys.

3) Robustness to some geometrical distortions: Although for the standard types of attacks (Table 2), the proposed method provides similar results as in [START_REF] Stanković | An Application of Multidimensional Time-Frequency Analysis as a base for the Unified Watermarking Approach[END_REF], it also provides certain robustness in the presence of geometrical attacks, such as image cropping and scaling with interpolation (the last two rows in Table 1), which is not achieved in the previous approach.

4) Security:

The watermark is shaped and added directly to the components in the space/spatialfrequency domain, and thus, it is hard to remove it without the key (which is assumed to be private). In other words, supposing that the image quality is important for the application, any attempt to remove the watermark will produce significant quality degradation. Furthermore, the method for pixel selection provides the candidates for watermark embedding that are equal from the perceptual point of view. Hence, generally, the watermark embedding can be done in any set of n points, where n<m (m is the number of candidates). These n pixels could be chosen by a random key, which provides the additional security of the procedure.

CONCLUSION

An image watermarking procedure for ownership protection that combines space/spatialfrequency analysis and Hermite projection method has been proposed. This approach offers several advantages regarding imperceptibility and robustness constraints, such as high global PSNR, higher local PSNR comparing to other space/spatial-frequency approaches and still low probability of error in the presence of attacks. Thus, it enhances the quality of watermarked image, providing better watermark imperceptibility. Each image region has its own space-spatialfrequency support which defines the watermark characteristics within the region, but the watermark detection is blind and does not need the information about support function.

Furthermore, although the Hermite expansion has been applied in both domains: spatial (for regions classification) and space/spatial-frequency domain (to create space/spatial-frequency mask), the complexity of the algorithm is not increased due to the existence of the fast Hermite projection method and fast recursive realization of Hermite functions.

Fig 1 .

 1 Fig 1. a) Block scheme of the watermarking procedure, b) Space/spatial-frequency

  denotes the STFT of the starting watermark created as a 2D pseudo-random sequence. In this way, the watermark will be separately adapted to the space/spatial-frequency content of each image region. In other words, each region has its own support function which is used to distribute watermark over the region in accordance with region components. This watermark design procedure should result in better watermark imperceptibility.

Fig 2 .

 2 Fig 2. a) An example of 2D space/spatial-frequency mask, b) watermark modelled using

  STFT of the original and watermarked image, respectively and α allows to tune the watermark strength. The value of α depends on the initial amplitude of pseudo-random sequence p. Since the amplitude of pseudo-random sequence influences the PSNR, alpha is set experimentally to achieve PSNR≈50 dB for a high quality of watermarked images.

Fig. 3 .Fig. 4 .

 34 Fig. 3. a) flat regions, b) moderately flat, c) moderately busy, d) busy

Fig 5 .Fig. 6 .

 56 Fig 5. An illustration of detector responses for 40 watermark (red) and 400 wrong trials

TABLE II .

 II COMPARISON WITH OTHER METHODS IN THE PRESENCE OF ATTACKS (IN TERMS OF PROBABILITY OF DETECTION ERROR)

	Type of Attack	Proposed Method in SSF domain	Previous Method in SSF domain	Standard DCT procedure	DWT domain procedure with GGF detector	DWT domain procedure with Cauchy detector
	No attacks	10 -15	10 -15	10 -8	10 -7	10 -9
	JPEG QF=80	10 -14	10 -13	10 -7	10 -6	10 -7
	JPEG QF=60	10 -12	10 -12	10 -5	10 -6	10 -6
	JPEG QF=40	10 -9	10 -7	10 -3	10 -5	10 -6
	Impulse noise 0.01	10 -9	10 -10	10 -6	10 -6	10 -6
	Gaussian noise 0.003	10 -9	10 -10	10 -6	10 -6	10 -5
	Darkening 40%	10 -14	10 -14	10 -7	10 -6	10 -7