[hdrcouchot.git] / stegoyousra.tex

La plupart des schémas de stéganographie sont conçus de sorte  à minimiser une 
fonction de distortion. Dans les exemples du chapitre précédent, 
ces fonctions de distortion sont construites dans l'obejctif de préserver 
les caractéristiques de l'images. 
On comprend aisément que dans des régions uniformes ou sur des bords clairement définis,
une modification même mineure de l'image est facilement détectable.
Au contraire dans les textures, le bruit ou les régions chaotiques 
sont  difficiles à modéliser. Les caractéristiques des images 
dont ces zones ont été modifiées sont ainsi similaires à celles
des images initiales.

Ces régions sont caractérisées par des courbes de niveau très perturbées.
Ce chapitre présente une nouvelle fonction de distortion pour la stéganographie 
qui est basée sur les dérivées du second ordre, l'outil mathématique usuel 
pour les courbes de niveau.

Pour peu qu'on sache définir une fonction $P$ 
qui associe à chaque pixel $(x,y)$ sa valeur $P(x,y)$,
les pixels tels que les dérivées secondes de $P$ ont des valeurs élevées 
sont des bon candidats pour contenir un bit du message.
Cependant, une telle fonction $P$ n'est connue que de manière discrète,
\textit{i.e.}, en un nombre fini de points. 
Les dérivées premières et secondes ne peuvent donc pas être évaluées mathématiquement.
Au mieux, on peut construire une fonction qui approxime ces $P$ sur cet ensemble 
de pixels. Ordonner alors les pixels selon la matrice hessienne 
(\textit{i.e.}, la matrice des dérivées secondes) n'est pas trivial puisque celle-ci 
contient de nombreuses valeurs pour un seul pixel donné. 

On verra dans ce chapitre comment des approximations des dérivées 
premières et secondes pour des images numriques (Section~\ref{sec:gradient}) on peu être
obtenues.
Deux propositions de dérivées secondes sont ensuite 
données et prouvées (Section~\ref{sec:second} et Section~\ref{sec:poly}).
Une adaptation d'une fonction de distortion existante est étudiée
en Section~\ref{sec:distortion} et des expériences sont présentées 
en Section~\ref{sec:experiments}.



\section{Des dérivées dans une image}\label{sec:gradient}

Cette section rappelle d'abord les liens entre lignes de niveau, gradient et 
matrice hessienne puis analyse ensuite leur construction à l'aide
de noyaux de la théorie du signal.


\subsection{Matrice hessienne}\label{sub:general}
On considère qu'une image peut être assimilée à une fonction de $\R^+\times \R^+$
dans $\R^+$ telle que la valeur $P(x,y)$ est associée à chaque  pixel de coordonnées $(x,y)$.
Les variations d'une telle fonction en $(x_0,y_0)$ peuvent être évaluées 
grace au gradient 
$\nabla{P}(x_0,y_0) = \left(\frac{\partial P}{\partial x}(x_0,y_0),\frac{\partial P}{\partial y}(x_0,y_0)\right).
$
Le vecteur gradient pointe dans la direction où la fonction a le plus fort acroissement.
Des pixels ayant des valeurs voisines sont sur des lignes de niveaux qui sont orthogonales
à ce vecteur.

Les variations du vecteur gradient s'expriment usuellement à l'aide de la matrice
hessienne $H$ des dérivées partielles de second ordre de $P$.
\[
H = \begin{bmatrix} 
\dfrac{\partial^2 P}{\partial x^2} & 
\dfrac{\partial^2 P}{\partial x \partial y} \\
\dfrac{\partial^2 P}{\partial y \partial x} & 
\dfrac{\partial^2 P}{\partial y^2} \\
\end{bmatrix}. 
\]

En un pixel $(x_0,y_0)$, plus les valeurs de cette matrice sont éloignées de zéro,
plus le gradient varie en ce point. Evaluer ce type de matrice est ainsi primordial
en stéganographie. Cependant cette tâche n'est pas aussi triviale qu'elle n'y 
paraît puisque les images naturelles ne sont pas  définies à l'aide 
de fonction différentiables de $\R^+\times \R^+$
dans $\R^+$. La suite montre comment obtenir des approximations de telles matrices. 

\subsection{Approches classiques pour évaluer le gradient dans des images}\label{sub:class:1}
Dans ce contexte, les approches les plus utilisées pour évaluer un gradient 
sont ``Sobel'', ``Prewitt'', ``Différence centrale'' et `` Difference intermediaire''.
Chacune de ces  approches applique un produit de convolution $*$ entre un noyau $K$
(rappelé dans le tableau~\ref{table:kernels:usual}) et une fenêtre $A$ de taille 
$3\times 3$. Le résultat 
 $A * K$ est une approximation du gradient horizontal
\textit{i.e.}, $\dfrac{\partial P}{\partial x}$.
Soit $K\rl$ le résultat de la rotation  d'un angle $\pi/2$ sur $K$.
La composante verticale du gradient,  $\dfrac{\partial P}{\partial y}$ est obtenue 
de manière similaire en évaluant $A * K\rl$. Lorsqu'on applique ceci sur toute 
la matrice image, on obtient peu ou prou une matrice de même taille pour chacune des 
dérivées partielles. 

Les deux éléments de la première ligne (respectivement de la seconde ligne) 
de la matrice hessienne
sont le résultat du calcul du gradient sur la matrice $\dfrac{\partial P}{\partial x}$
(resp. sur la matrice  $\dfrac{\partial P}{\partial y}$).


\begin{table}[ht]
\caption{Noyaux usuels pour évaluer des gradients d'images\label{table:kernels:usual}
}
\centering
\scriptsize
\begin{tabular}{|c|c|c|}
    \hline
    Nom&   Sobel & Prewitt \\
    \hline
    Noyau & $\textit{Ks}= \begin{bmatrix} -1 & 0 & +1 \\ -2 & 0 & +2 \\ -1 & 0 & +1 \end{bmatrix} $ &
    $\textit{Kp}= \begin{bmatrix} -1 & 0 & +1 \\ -1 & 0 & +1 \\ -1 & 0 & +1 \end{bmatrix} $\\
    \hline\hline
    Nom & Difference  & Différence  \\
            & centrale & Intermediaire \\ 
    \hline
    Noyau & $\textit{Kc}= \begin{bmatrix} 0&0&0 \\ -\dfrac{1}{2} & 0 & +\dfrac{1}{2} \\ 0&0&0 \end{bmatrix} $ &
    $\textit{Ki}= \begin{bmatrix} 0&0&0 \\ 0 & -1 & 1 \\ 0&0&0 \end{bmatrix} $ \\
    \hline
\end{tabular}
\end{table}




\subsection{Matrices hessiennes induites  par des approches 
de gradient d'images}\label{sub:class:2}
Il est connu que  
$\dfrac{\partial^2 P}{\partial x \partial y} $ est égal à 
$\dfrac{\partial^2 P}{\partial y \partial x}$ si 
les méthode qui calculent le grandient et le gradient du gradient (la matrice hessienne)
sont les mêmes.

For instance, in the Sobel approach, 
it is easy to verify that the calculus of 
$\dfrac{\partial^2 P}{\partial x \partial y}$ 
and of 
$\dfrac{\partial^2 P}{\partial y \partial x}$ 
are both the result 
of a convolution product with the Kernel 
$\textit{Ks}''_{xy}$ given in Table~\ref{table:hessian:usual}.
This one summarizes kernels
$K_{x^2}''$ and 
$K_{xy}''$  
that allow to respectively compute  
$\dfrac{\partial^2 P}{\partial x^2}$ and
$\dfrac{\partial^2 P}{\partial x \partial y}$ with a convolution product
for each of the usual image gradient operator.

\begin{table}[ht]
\caption{Kernels of second order gradient operators\label{table:hessian:usual}
}\centering
\tiny
\begin{tabular}{|c|c|}
    \hline
    Sobel & Prewitt \\
    \hline
    $
    \textit{Ks}_{x^2}''= 
    \begin{bmatrix} 
    1 &  0  & -2  & 0 & 1 \\ 
    4 &  0  & -8  & 0 & 4 \\ 
    6 &  0  & -12 & 0 & 6 \\
    4 &  0  & -8  & 0 & 4 \\ 
    1 &  0  & -2  & 0 & 1 
    \end{bmatrix} 
    $
    & 
    $
    \textit{Kp}_{x^2}''= 
    \begin{bmatrix} 
    1 &  0  & -2  & 0 & 1 \\ 
    2 &  0  & -4  & 0 & 2 \\ 
    3 &  0  & -6 & 0 & 3 \\
    2 &  0  & -4  & 0 & 2 \\ 
    1 &  0  & -2  & 0 & 1 
    \end{bmatrix} 
    $
    \\
    \hline
    $
    \textit{Ks}_{xy}''=
    \begin{bmatrix} 
    -1 & -2 & 0 & 2 & 1 \\ 
    -2 & -4 & 0 & 4 & 2 \\ 
    0  & 0  & 0 & 0 & 0 \\
    2 & 4 & 0 & -4 & -2 \\ 
    1 & 2 & 0 & -2 & -1 
    \end{bmatrix} 
    $
    &
    $
    \textit{Kp}_{xy}''=
    \begin{bmatrix} 
    -1 & -1 & 0 & 1 & 1 \\ 
    -1 & -1 & 0 & 1 & 1 \\ 
    0  & 0  & 0 & 0 & 0 \\
    1 & 1 & 0 & -1 & -1 \\ 
    1 & 1 & 0 & -1 & -1 
    \end{bmatrix} 
    $ \\
    
    \hline\hline
    Central & Intermediate \\
    Difference &Difference \\ 
    \hline
    $
    \textit{Kc}_{x^2}''= 
    \begin{bmatrix} 
    0 &  0  & 0  & 0 & 0 \\ 
    0 &  0  & 0  & 0 & 0 \\ 
    \dfrac{1}{4} &  0  & -\dfrac{1}{2} & 0 & \dfrac{1}{4} \\
    0 &  0  & 0  & 0 & 0 \\ 
    0 &  0  & 0  & 0 & 0  
    \end{bmatrix} 
    $
    & 
    $
    \textit{Ki}_{x^2}''= 
    \begin{bmatrix} 
    0 &  0  & 0  & 0 & 0 \\ 
    0 &  0  & 0  & 0 & 0 \\ 
    0 &  0  & 1  & -2 & 1 \\
    0 &  0  & 0  & 0 & 0 \\ 
    0 &  0  & 0  & 0 & 0  
    \end{bmatrix} 
    $
    \\
    \hline
    $
    \textit{Kc}_{xy}''= 
    \begin{bmatrix} 
     -\dfrac{1}{4}  & 0 & \dfrac{1}{4}\\ 
     0  & 0 & 0 \\ 
     \dfrac{1}{4} & 0 &  -\dfrac{1}{4} 
    \end{bmatrix} 
    $
    & 
    $
    \textit{Ki}_{xy}''= 
    \begin{bmatrix} 
     0  & -1 & 1 \\ 
     0  & 1 & -1 \\ 
     0  &  0 & 0 
    \end{bmatrix} 
    $\\
    \hline
\end{tabular}

\end{table}

The Sobel kernel $\textit{Ks}_{x^2}''$ allows to detect whether the central 
pixel belongs to a ``vertical'' edge, even if this one is noisy, by considering its 
vertical neighbours. The introduction of these vertical neighbours in this kernel is meaningful 
in the context of finding edges, but not very accurate when the objective is to 
precisely find the level curves of the image. 
Moreover, all the pixels that are in the  second and the fourth column in this kernel
are ignored.
The Prewitt Kernel has similar drawbacks in this context.


The Central Difference kernel $\textit{Kc}_{x^2}''$ is not influenced by 
the vertical neighbours of the central pixel and is thus more accurate here.
However, the kernel $\textit{Kc}_{xy}''$ again looses the values of the pixels that 
are vertically and diagonally aligned with the central one.

Finally, the Intermediate Difference kernel $\textit{Ki}_{x^2}''$ shifts 
to the left the value of horizontal variations of $\dfrac{\partial P}{\partial x}$:
the central pixel $(0,0)$ exactly receives the value
$\dfrac{P(0,2)-P(0,1)}{1} - \dfrac{P(0,1)-P(0,0)}{1}$,
which is an approximation of 
$\dfrac{\partial P}{\partial x}(0,1)$ and not of 
$\dfrac{\partial P}{\partial x}(0,0)$.
Furthermore the Intermediate Difference kernel $\textit{Ki}_{xy}''$ only deals 
with pixels in the upper right corner, loosing all the other information.

Due to these drawbacks, we are then left to produce another approach to find the level curves with strong accuracy.

\section{Second Order Kernels for Accurate Level Curves}\label{sec:second}

This step aims at finding 
accurate level curve variations in an image. 
We do not restrict the kernel to have a fixed size (\textit{e.g.}, $3\times3$ or $5 \times 5$ as in the 
aforementioned schemes).
This step is thus defined with kernels of size 
$(2n+1)\times (2n+1)$, $n \in \{1,2,\dots,N\}$, where 
$N$ is a parameter of the approach.

The horizontal gradient variations are thus captured thanks to $(2n+1)\times (2n+1)$ square kernels
\begin{small}
\[
\arraycolsep=1.4pt
\def\arraystretch{1.4}
\def\arraystretch{1.4}
    \textit{Ky}_{x^2}''=
    \left(
    \begin{array}{ccccccccc}
         0           & & & & \dots& & & & 0  \\
         \vdots      & & & & & & & & \vdots \\
         0           & & & & \dots & & & & 0  \\
         \dfrac{1}{2n}& 0 & \dots & 0  & -\dfrac{2}{2n} & 0  & \dots & 0& \dfrac{1}{2n} \\
         0           & & & & \dots& & & & 0  \\
        \vdots      & & & & & & & & \vdots \\
         0           & & & & \dots & & & & 0  
    \end{array}
    \right)
\]
\end{small}

When the convolution product is applied on a $(2n+1)\times(2n+1)$ window,
the result is 
$\dfrac{1}{2}\left(\dfrac{P(0,n)-P(0,0)}{n} - \dfrac{P(0,0)-P(0,-n)}{n}\right)$, which is indeed 
the variation between the gradient around the central pixel. 
This proves that this calculus is a correct approximation of 
$\dfrac{\partial^2 P}{\partial x^2}$.

When $n$ is 1, this kernel is a centered version of the horizontal Intermediate Difference kernel $\textit{Ki}_{x^2}''$ modulo a multiplication by $1/2$. When $n$ is 2, this kernel is equal to $\textit{Kc}_{x^2}''$.

The vertical gradient variations are again obtained by applying 
a $\pi/2$ rotation to each horizontal kernel $\textit{Ky}_{x^2}''$.

The diagonal gradient variations are obtained thanks to the $(2n+1)\times (2n+1)$ square kernels 
$\textit{Ky}_{xy}''$ defined by

\begin{small}
\[
\arraycolsep=1.4pt
\def\arraystretch{1.4}
\textit{Ky}_{xy}'' = \dfrac{1}{4}
\left(
    \begin{array}{ccccccccc}
     \frac{1}{n^2}& \dots & \frac{1}{2n} & \frac{1}{n} 
     & 0 &
     -\frac{1}{n}&-\frac{1}{2n} & \dots &  -\frac{1}{n^2} 
     \\
     \vdots & 0     &    &
     & \dots &
      &  &  0 & \vdots 
     \\
     \frac{1}{2n} & 0 &        &
     & \dots &
        &  & 0& -\frac{1}{2n} 
     \\
     \frac{1}{n} & 0 &    &    
     & \dots &
      &  & 0 &  -\frac{1}{n}
     \\
     0      &      & & & \dots& & & & 0  \\
     -\frac{1}{n} & 0 &        &
     & \dots &
        &  &0 & \frac{1}{n}
     \\
     -\frac{1}{2n} & 0 &       &
     & \dots &
      &  & 0 &  \frac{1}{2n} 
     \\
      \vdots & 0 &    &    
     & \dots &
        &  & 0& \vdots 
     \\
     -\frac{1}{n^2}&  \dots & -\frac{1}{2n} & -\frac{1}{n} 
     & 0 &
     \frac{1}{n}& \frac{1}{2n} & \dots  & \frac{1}{n^2}
    \end{array}
    \right).
\]
\end{small}

When $n$ is 1, $\textit{Ky}_{xy}''$ is equal to the kernel
$\textit{Kc}_{xy}''$, and %.
%
%When $n$ is 1, 
the average vertical variations of the horizontal variations are 
\[
\begin{array}{l}
\dfrac{1}{4}
\left[ 
\left((P(0,1)-P(0,0))-(P(1,1)-P(1,0))\right)+ \right.\\
\quad \left((P(-1,1)-P(-1,0))-(P(0,1)-P(0,0))\right)+\\
\quad \left((P(0,0)-P(0,-1))-(P(1,0)-P(1,-1))\right)+\\
\quad \left. \left((P(-1,0)-P(-1,-1))-(P(0,0)-P(0,-1))\right)
\right]  \\
=\dfrac{1}{4}
\left[ P(1,-1) -P(1,1) - P(-1,-1) + P(-1,1)\right].
\end{array}
\]
which is $\textit{Ky}_{xy}''$.

Let us now consider any number $n$, $1 \le n \le N$.
Let us first investigate the vertical variations related to 
the horizontal vector $\vec{P_{0,0}P_{0,1}}$ 
(respectively  $\vec{P_{0,-1}P_{0,0}}$)
of length 1 that starts from (resp. that points to) $(0,0)$. 
As with the case $n=1$, there are 2 new vectors of 
length 1, namely 
$\vec{P_{n,0}P_{n,1}}$ and $\vec{P_{-n,0}P_{-n,1}}$ 
(resp. 
$\vec{P_{n,-1}P_{n,0}}$, and $\vec{P_{-n,-1}P_{-n,0}}$)
that are vertically aligned with $\vec{P_{0,0}P_{0,1}}$ 
(resp. with $\vec{P_{0,-1}P_{0,0}}$).

The vertical variation is now equal to $n$. Following the case where $n$ is 1 to compute the average variation, 
the coefficients of the first and last line around the central 
vertical line are thus from left to right:
$\dfrac{1}{4n}$,
$\dfrac{-1}{4n}$,
$\dfrac{-1}{4n}$, and
$\dfrac{1}{4n}$.

Cases are similar with vectors $\vec{P_{0,0}P_{0,1}}$, \ldots  
$\vec{P_{0,0}P_{0,n}}$ which respectively lead to coefficients 
$-\dfrac{1}{4 \times 2n}$, \ldots, 
$-\dfrac{1}{4 \times n.n}$, and the proof is omitted.
Finally, let us consider the vector $\vec{P_{0,0}P_{0,1}}$
and its vertical variations when $\delta y$ is $n-1$.
As in the case where $n=1$, we thus obtain the coefficients 
$\dfrac{1}{4 \times (n-1)n}$ and 
$-\dfrac{1}{4 \times (n-1)n}$ 
 (resp. $-\dfrac{1}{4 \times (n-1)n}$ and 
$\dfrac{1}{4 \times (n-1)n}$)
in the second line (resp. in the 
penultimate line) since the vector has length $n$
and $\delta y$ is $n-1$.
Coefficient in the other lines are similarly obtained and the proof is thus omitted.

We are then left to compute an approximation of the partial second order derivatives 
$\dfrac{\partial^2 P}{\partial x^2}$,  $\dfrac{\partial^2 P}{\partial y^2}$, and $\dfrac{\partial^2 P}{\partial x \partial y}$
with the kernels, 
$\textit{Ky}_{x^2}''$, $\textit{Ky}_{y^2}''$, and $\textit{Ky}_{xy}''$ respectively.
However, the size of each of these kernels is varying from $3\times3$ to $(2N+1)\times (2N+1)$.
Let us explain the approach on the former partial derivative.
The other can be immediately deduced.

Since the objective is to detect large variations, the second order derivative is approximated as 
the maximum of the approximations. More formally, 
let $n$, $1 \le n \le N$, be an integer number and 
$\dfrac{\partial^2 P}{\partial x^2}_n$ be the result of applying the Kernel 
$\textit{Ky}_{x^2}''$ of size $(2n+1)\times (2n+1)$. The derivative 
$\dfrac{\partial^2 P}{\partial x^2}$ is defined by

\begin{equation}
\dfrac{\partial^2 P}{\partial x^2}
= \max \left\{ 
\abs{\dfrac{\partial^2 P}{\partial x^2}_1}, \dots, \abs{\dfrac{\partial^2 P}{\partial x^2}_N}
\right \}. 
\label{eq:d2p_dx2}
\end{equation}

The same iterative approach is applied to compute approximations of
$\dfrac{\partial^2 P}{\partial y \partial x}$ 
and of 
$\dfrac{\partial^2 P}{\partial y^2}$.
Next section studies the suitability of approximating second order derivatives
when considering an image as a polynomial.


\section{Polynomial Interpolation of Images for Hessian Matrix Computation}\label{sec:poly}
Let $P(x,y)$ be the discrete value of the pixel $(x,y)$ in the image.
Let $n$, $1 \le n \le N$, be an integer such that the objective is to find a polynomial interpolation 
on the $(2n+1)\times(2n+1)$ window where the central pixel has index $(0,0)$.  
There exists an unique polynomial $L : \R\times \R \to \R$ of degree $(2n+1)\times(2n+1)$ defined  
such that $L(x,y)=P(x,y)$ for each pixel $(x,y)$ in this window.
Such a polynomial is defined by 
\begin{equation}
\begin{array}{l}
L(x,y) =  
\sum_{i=-n}^{n}
\sum_{j=-n}^{n} \\
\quad P(i,j)
\left( 
\prod_{\stackrel{-n\leq j'\leq n}{j'\neq j}}
\frac{x-j'}{i-j'}
\right)
\left( 
\prod_{\stackrel{-n\leq i'\leq n}{i'\neq i}}
\frac{x-i'}{i-i'}
\right)
\end{array}
\end{equation}

It is not hard to prove that the first order horizontal derivative of the polynomial $L(x,y)$ 
is
\begin{equation}
\begin{array}{l}
\dfrac{\partial L}{\partial x} =  
\sum_{i=-n}^{n}
\sum_{j=-n}^{n} 
P(i,j)
\left( 
\prod_{\stackrel{-n\leq j'\leq n}{j'\neq j}}
\frac{y-j'}{j-j'}
\right)\\
\quad
\left( 
\sum_{\stackrel{-n\leq i'\leq n}{i'\neq i}}
\frac{1}{i-i'}
\prod_{\stackrel{-n\leq i''\leq n}{i''\neq i,i'}}
\frac{x-i''}{i-i''}
\right)
\end{array}
\end{equation}
\noindent and thus to deduce that the
second order ones are

\begin{equation}
\begin{array}{l}
\dfrac{\partial^2 L}{\partial x^2} =  
\sum_{i=-n}^{n}
\sum_{j=-n}^{n} 
P(i,j)
\left( 
\prod_{\stackrel{-n\leq j'\leq n}{j'\neq j}}
\frac{y-j'}{j-j'}
\right)\\
\quad
\left( 
\sum_{\stackrel{-n\leq i'\leq n}{i'\neq i}}
\frac{1}{i-i'}
\sum_{\stackrel{-n\leq i''\leq n}{i''\neq i,i'}}
\frac{1}{i-i''}
\prod_{\stackrel{-n\leq i'''\leq n}{i'''\neq i,i',i''}}
\frac{x-i'''}{i-i'''}
\right)
\end{array}
\label{eq:deriv:poly:x2}
\end{equation}

\begin{equation}
\begin{array}{l}
\dfrac{\partial^2 L}{\partial y \partial x} =  
\sum_{i=-n}^{n}
P(i,j) \\
\quad
\left( 
\sum_{\stackrel{-n\leq j'\leq n}{j'\neq j}}
\frac{1}{j-j'}
\prod_{\stackrel{-n\leq j''\leq n}{j''\neq j, j'}}
\frac{y-j''}{j-j''}
\right)\\
\quad
\left( 
\sum_{\stackrel{-n\leq i'\leq n}{i'\neq i}}
\frac{1}{i-i'}
\prod_{\stackrel{-n\leq i''\leq n}{i''\neq i, i'}}
\frac{x-i''}{i-i''}
\right)
\end{array}
\label{eq:deriv:poly:yx}
\end{equation}

These second order derivatives  are computed for each moving 
window and are associated to the central pixel, \textit{i.e.}, to the pixel $(0,0)$ inside this one.

Let us first simplify $\dfrac{\partial^2 L}{\partial x^2}$ when $(x,y)=(0,0)$
defined in Equation~(\ref{eq:deriv:poly:x2}). If $j$ is not null, the index $j'$
is going to be null and the product  
$\left( 
\prod_{\stackrel{-n\leq j'\leq n}{j'\neq j}}
\frac{-j'}{j-j'}
\right)$ is null too. 
In this equation, we thus only consider $j=0$.
It is obvious that the product indexed with $j'$ is thus equal to 1.
This equation can thus be simplified in:

\begin{equation}
\begin{array}{l}
\dfrac{\partial^2 L}{\partial x^2} =  
\sum_{i=-n}^{n}
P(i,0)\\
\quad
\left( 
\sum_{\stackrel{-n\leq i'\leq n}{i'\neq i}}
\frac{1}{i-i'}
\sum_{\stackrel{-n\leq i''\leq n}{i''\neq i,i'}}
\frac{1}{i-i''}
\prod_{\stackrel{-n\leq i'''\leq n}{i'''\neq i,i',i''}}
\frac{i'''}{i'''-i}
\right)
\end{array}
\label{eq:deriv:poly:x2:simpl}
\end{equation}


and then in:

\begin{equation}
\begin{array}{l}
\dfrac{\partial^2 L}{\partial x^2} =  
\sum_{i=-n}^{n}
P(i,0)\\
\quad
\left( 
\sum_{\stackrel{-n\leq i' < i'' \le n}{i',i''\neq i}}
\frac{2}{(i-i')(i-i'')}
\prod_{\stackrel{-n\leq i'''\leq n}{i'''\neq i,i',i''}}
\frac{i'''}{i'''-i}
\right).
\end{array}
\label{eq:deriv:poly:x2:simpl:2}
\end{equation}

From this equation, the kernel allowing to evaluate horizontal 
second order derivatives 
can be computed for any $n$.
It is further denoted as $Ko''_{x^2}$. 
Instances of such matrix when $n=2$, $3$, and $4$
are given in Table~\ref{table:sod:hori:poly}.

\begin{table}[ht]
    \caption{Kernels $Ko''_{x^2}$ 
    for second order horizontal derivatives induced 
    by polynomial interpolation}
    \centering
    \scriptsize
\def\arraystretch{1.4}
    \begin{tabular}{|c|c|}
         \hline
         $n$ & $Ko''_{x^2}$ \\
         \hline
         $2$ & $\left[\dfrac{-1}{12}, \dfrac{4}{3} , \dfrac{-5}{2}, \dfrac{4}{3} \dfrac{-1}{12}\right]$ \\
         \hline
         $3$ & $\left[\dfrac{1}{90}, \dfrac{-3}{20}, \dfrac{3}{2}, \dfrac{-49}{18}, \dfrac{3}{2}, \dfrac{-3}{20}, \dfrac{1}{90}\right]$ \\
         \hline
         $4$ & $\left[\dfrac{-1}{560}, \dfrac{8}{315}, \dfrac{-1}{5}, \dfrac{8}{5}, \dfrac{-205}{72}, \dfrac{8}{5}, \dfrac{-1}{5}, \dfrac{8}{315}, \dfrac{-1}{560}\right]$\\
         \hline
         %$5$ & $\left[\dfrac{1}{3150}, \dfrac{-5}{1008}, \dfrac{5}{126}, %\dfrac{-5}{21}, \dfrac{5}{3}, \dfrac{-5269}{1800}, \dfrac{5}{3}, %\dfrac{-5}{21}, \dfrac{5}{126}, \dfrac{-5}{1008}, \dfrac{1}{3150}\right]$\\
         %\hline
    \end{tabular}

    \label{table:sod:hori:poly}
\end{table}


\begin{table}[ht]
\caption{Kernels for second order diagonal derivatives induced 
    by polynomial interpolation \label{table:sod:diag:poly}
}
\centering
\scriptsize
\def\arraystretch{1.5}
\begin{tabular}{|c|c|}
    \hline
    $n$ & $Ko''_{xy}$\\
    \hline
    2 & $
\begin{bmatrix} 
\dfrac{1}{4} & 0 &  \dfrac{-1}{4}\\
 0 & 0 &0\\
\dfrac{-1}{4} & 0 &  \dfrac{1}{4}\\
\end{bmatrix}
$
 \\   
 \hline
   3 &  
   $\begin{bmatrix} 
\dfrac{1}{144} & \dfrac{-1}{18} & 0 & \dfrac{1}{18} & \dfrac{-1}{144}\\
\dfrac{-1}{18} & \dfrac{4}{9} & 0 & \dfrac{-4}{9} & \dfrac{1}{18}\\
0 & 0 & 0 & 0 &0\\
\dfrac{1}{18} & \dfrac{-4}{9} & 0 & \dfrac{4}{9} & \dfrac{-1}{18}\\
\dfrac{-1}{144} & \dfrac{1}{18} & 0 & \dfrac{-1}{18} & \dfrac{1}{144}
\end{bmatrix}
$\\
  \hline
\end{tabular}
\end{table}


From Equation~(\ref{eq:deriv:poly:yx}), kernels allowing to evaluate diagonal 
second order derivatives (\textit{i.e.}, 
$\dfrac{\partial^2 L}{\partial y \partial x}$)  are computed.
They are denoted as $Ko''_{xy}$.
Table~\ref{table:sod:diag:poly} gives two examples of them when $n=1$ and $n=2$.
Notice that for $n=1$, the kernel $Ko''_{xy}$ is equal to $Kc''_{xy}$. 






\section{Distortion Cost}\label{sec:distortion}

The distortion function has to associate to each pixel $(i,j)$ 
the cost $\rho_{ij}$ of its modification
by $\pm 1$.

The objective is to map a small value to a pixel when all its second order derivatives 
are high and a large value otherwise. 
%\YF{The WOW  use the  wavelet relative distortion (UNIWARD) as a function to force the embedding operation to the noisy areas and avert the smooth areas.}
In WOW and UNIWARD  the distortion function is based on the H\"older norm with
\[
\rho_{ij}^w = 
\left(
\abs{\xi_{ij}^h}^{p} +
\abs{\xi_{ij}^v}^{p} +
\abs{\xi_{ij}^d}^{p} 
\right)^{-\frac{1}{p}}
\]
where $p$ is a negative number and 
$\xi_{ij}^h$ (resp. $\xi_{ij}^v$ and $\xi_{ij}^d$)
represents the horizontal (resp. vertical and diagonal) suitability.
A small suitability in one direction means an inaccurate position to embed a message.

We propose here to adapt such a distortion cost as follows:
\[
\rho_{ij} = 
\left(
\abs{\dfrac{\partial^2 P}{\partial x^2}(i,j)} +
\abs{\dfrac{\partial^2 P}{\partial y^2}(i,j)} +
\abs{\dfrac{\partial^2 P}{\partial y \partial x}(i,j)}
\right)^{-\frac{1}{p}}
\]
It is not hard to check that such a function has large 
value when at least one of its derivatives is null. Otherwise,
the larger the derivatives are, the smaller the returned value
is.




%Practically, it is computed as a double convolution product with Daubechies Wavelet Kernels. 

%It is not hard to see that this approach 
%does not consider the range of the values. For %instance, if one direction only has large suitability
%absolute values, all the elements $\abs{\xi_{ij}}^{p}$ are so small, and 
%suitable pixels cannot be distinguished from inaccurate ones. The distortion function has to 
%provide a way to discriminate pixels even when the distribution has small deviation

%Because of each pixel value ranges in $[0,255]$, each first order derivative belongs to
%$[-255,255]$ and thus each second order derivative has a real value in 
%$[-255,255]$.  
%Due to the definition of the approximation of $\dfrac{\partial^2 P}{\partial x^2}$, 
%this one is a positive or null real number lower than 255. 

%The higher the second order derivatives are, 
%the 

%To enlarge the discrimination power, 
%data are firstly linearly mapped to $[0,1]$. 
%The tuple 
%$(\dfrac{\partial^2 P}{\partial x^2} , \dfrac{\partial^2 P}{\partial y^2}, \dfrac{\partial^2 %P}{\partial x \partial y})$ belongs to $[0,1]\times[0,1]\times[0,1] = [0,1]^3$.

%Basically, for each dimension the idea is to select a subset of pixels of size $k$ and next to retain those which have been selected three times, \textit{i.e.}, those with large values 
%in any element of 
%$(\dfrac{\partial^2 P}{\partial x^2} , \dfrac{\partial^2 P}{\partial y^2}, \dfrac{\partial^2 %P}{\partial x \partial y})$.
%Since the dimension is $3$, the size $k$ can be defined by
%\[
%k = W \times H \times \alpha ^{\dfrac{1}{3}}
%\]
%where $W$ and $H$ respectively stand for the width and the length of the image (expressed in pixels)
%and $\alpha$ is the payload. For instance if $\alpha$ is 0.1, 
%about $0.46\%$ of the image pixels with higher values in $\dfrac{\partial^2 P}{\partial x^2}$ should be selected, and similarly with $\dfrac{\partial^2 P}{\partial y^2}$ and $\dfrac{\partial^2  P}{\partial x \partial y}$.


% Dire comment on les ordonne, les autres et le sens

\begin{figure*}
    \centering
    \begin{tabular}{|c|c|c|}
    \hline
    Scheme & Stego. content & Changes with cover \\
    \hline
    $Ky$ based approach &\includegraphics[scale=0.20]{images/38_dp}&
    \includegraphics[scale=0.20]{images/38_dp_diff} \\
    \hline
    $Ko$ based approach & \includegraphics[scale=0.20]{images/38_poly} &
    \includegraphics[scale=0.20]{images/38_poly_diff} \\
    \hline
    \end{tabular}
    \caption{Embedding changes instance with payload  $\alpha = 0.4$}
    \label{fig:oneimage}
\end{figure*}


\section{Experiments}\label{sec:experiments}

First of all, the whole steganographic approach code is available online\footnote{\url{https://github.com/stego-content/SOS}}.

Figure~\ref{fig:oneimage} presents the results of embedding data in a cover image
from the BOSS contest database~\cite{Boss10} with respect to the two second order derivative schemes presented in this work. 
The $Ky$ based approach (resp. the $Ko$ based one)
corresponds to the scheme detailed in Section~\ref{sec:second} 
(resp. in Section~\ref{sec:poly}).
The payload $\alpha$ is set to 0.4 and kernels are computed with $N=4$.    
The central column outputs the embedding result whereas the right one displays 
differences between the cover image and the stego one.
It can be observed that pixels in smooth area (the sky, the external access steps) and pixels in clean edges (the columns, 
the step borders) are not modified by the approach.
On the contrary, an unpredictable area (a monument for example) concentrates pixel changes.


\subsection{Choice of parameters}

The two methods proposed in Section~\ref{sec:second} and  
in Section~\ref{sec:poly} are based on kernels of size up to 
$(2N+1)\times(2N+1)$. This section aims at finding the 
value of the $N$ parameter that maximizes the security level.
For each approach, we have built 1,000 stego images with 
$N=2$, $4$, $6$, $8$, $10$, $12$, and $14$ where the covers belong 
to the BOSS contest database. 
This set contains 10,000 grayscale $512\times 512$ images in a RAW format. 
The security of the approach has been evaluated thanks to the Ensemble Classifier~\cite{DBLP:journals/tifs/KodovskyFH12} based steganalyser,
which is considered as a state of the art steganalyser tool.
This steganalysis process embeds the rich model (SRM) features~\cite{DBLP:journals/tifs/FridrichK12}
of size 34,671. 
For a payload $\alpha$, either equal to $0.1$ or to $0.4$, average testing errors (expressed in percentages) have been studied and are summarized in
Table~\ref{table:choice:parameter}.

\begin{table}[ht]
\caption{Average Testing Errors with respect to the the Kernel Size}
\begin{scriptsize}
\centering
\setlength{\tabcolsep}{3pt} 
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
\cline{2-9}
\multicolumn{1}{c|}{} & \multirow{2}{*}{$\alpha$} & \multicolumn{7}{c|}{$N$} \\
\cline{3-9}
\multicolumn{1}{c|}{}&  & $2$ & $4$&  $6$&  $8$& $10$& $12$ & $14$ \\
\hline{}
Average testing  
& \textit{0.1} & 39& 40.2&  39.7&  39.8& 40.1& $39.9$&  $39.8$  \\
\cline{2-9}
error for Kernel $K_y$ & \textit{0.4}& 15& 18.8& 19.1&  19.0& 18.6& 18.7 & 18.7 \\
\hline
Average testing & \textit{0.1} & 35.2 & 36.6&  36.7&  36.6& 37.1& 37.2 & 37.2 \\
\cline{2-9}
error for Kernel $K_o$ & \textit{0.4}  & 5.2 & 6.8& 7.5 & 7.9 & 8.1 & 8.2 & 7.6 \\
\hline
\end{tabular}
\label{table:choice:parameter}
\end{scriptsize}
\end{table}

Thanks to these experiments, we observe that the size $N=4$ (respectively $N=12$) obtains sufficiently  large average testing errors for the $Ky$ based approach 
(resp. for the $Ko$ based one). In what follows, these values are retained for these two methods.

\subsection{Security Evaluation}
As in the previous section, the BOSS contest database has been retained.
To achieve a complete comparison with other steganographic tools,
the whole database of 10,000 images has been used.
Ensemble Classifier with SRM features is again used to evaluate the security of the approach.

We have chosen 4 different payloads, 0.1, 0.2, 0.3, and 0.4, as in many steganographic evaluations.
Three values are systematically given for each experiment:
the area under the ROC curve (AUC),
the average testing error (ATE),
and the OOB error (OOB).

All the results are summarized in Table~\ref{table:experiments:summary}.
%and  in Figure~\ref{fig1}, Figure~\ref{fig2}, and Figure~\ref{fig3}.
Let us analyse these experimental results. 
The security approach is often lower than those observed with state of the art tools:
for instance with payload $\alpha=0.1$, the most secure approach is WOW 
with an average testing error equal to 0.43 whereas our approach reaches 0.38.
However these results are  promising and for two reasons.
First, our approaches give more resistance towards Ensemble Classifier (contrary to HUGO)
for large payloads.
Secondly, without any optimisation, our approach is not so far from state of the art steganographic tools.
Finally, we explain the lack of security of the $Ko$ based approach with large payloads as follows:
second order derivatives are indeed directly extracted from polynomial interpolation.
This easy construction however induces large variations between the polynomial $L$ and 
the pixel function $P$.

\begin{table}
\caption{Summary of experiments}\label{table:experiments:summary}
\begin{small}
\centering
\begin{tabular}{|l|l|l|l|l|}
\hline

 & Payload & AUC & ATE   &  OOB \\ \hline
{WOW}   
            & 0.1 & 0.6501 & 0.4304 & 0.3974\\
            & 0.2 & 0.7583 & 0.3613 & 0.3169\\
            & 0.3 & 0.8355 & 0.2982 & 0.2488\\
            & 0.4 & 0.8876 & 0.2449 & 0.1978\\ 
                                       \hline

{SUNIWARD}  & 0.1 & 0.6542 & 0.4212 & 0.3972\\
            & 0.2 & 0.7607 & 0.3493 & 0.3170\\
            & 0.3 & 0.8390 & 0.2863 & 0.2511\\
            & 0.4 & 0.8916 & 0.2319 & 0.1977\\ 
 \hline
 {MVG}      & 0.1 & 0.6340 & 0.4310 &0.4124 \\
            & 0.2 & 0.7271 & 0.3726 &0.3399 \\
            & 0.3 & 0.7962 & 0.3185& 0.2858\\
            & 0.4 & 0.8486& 0.2719& 0.2353 \\ 
 \hline
  {HUGO}    & 0.1 & 0.6967 & 0.3982 & 0.3626 \\
            & 0.2 & 0.8012 & 0.3197 & 0.2847 \\
            & 0.3 & 0.8720 & 0.2557 & 0.2212 \\
            & 0.4 & 0.9517 & 0.1472 & 0.1230 \\ 
 \hline
{$Ky$ based approach} 
            & 0.1 & 0.7378 & 0.3768 & 0.3306 \\
            & 0.2 & 0.8568 & 0.2839 & 0.2408 \\
            & 0.3 & 0.9176 & 0.2156 & 0.1710 \\
            & 0.4 & 0.9473 & 0.1638 & 0.1324\\
 \hline
{$Ko$ based approach} 
            & 0.1 & 0.6831 & 0.3696  & 0.3450 \\
            & 0.2 & 0.8524 & 0.1302  & 0.2408 \\
            & 0.3 & 0.9132 & 0.1023  & 0.1045 \\
            & 0.4 & 0.9890 & 0.0880  & 0.0570 \\
                              \hline 
                              
\end{tabular}
\end{small}
\end{table}


%\begin{figure}[ht]
%\centering
%\includegraphics[width=9cm]{22} 
%\caption{Average Testing Error(ATE) } 
%\label{fig1}
%\end{figure}

%\begin{figure}[ht]
%\centering
%\includegraphics[width=9cm]{11} 
%\caption{Out of Bag($E_{\textit{OOB}}$)} 
%\label{fig2}
%\end{figure}

%\begin{figure}[ht]
%\centering
%\includegraphics[width=9cm]{33} 
%\caption{Area Under Curve(AUC)} %\label{fig3}
%\end{figure}


\section{Conclusion}

The first contribution of this paper is to propose of a distortion
function which is based on second order derivatives. These
partial derivatives allow to accurately compute 
the level curves and thus to look favorably on pixels
without clean level curves. 
Two approaches to build these derivatives have been proposed.
The first one is based on revisiting kernels usually embedded 
in edge detection algorithms. 
The second one is based on the polynomial approximation
of the bitmap image.
These two methods have been completely implemented.
The first experiments have shown that the security level 
is slightly inferior the one of the most stringent approaches. These first promising results encourage us to deeply investigate this research direction. 

Future works aiming at improving the security of this proposal are planned as follows. The authors want first to focus on other approaches to provide second order derivatives with larger discrimination power.
Then, the objective will be to deeply investigate whether the H\"older norm is optimal when the objective is to avoid null second order derivatives, and to give priority to the largest second order values.