Installing toolboxes and setting up the path.

You need to unzip these toolboxes in your working directory, so that you have toolbox_signal and toolbox_general in your directory.

For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.

Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.

Execute this line only if you are using Matlab.

getd = @(p)path(p,path); % scilab users must *not* execute this

Then you can add the toolboxes to the path.

getd('toolbox_signal/');
getd('toolbox_general/');

Digital Image Watermarking

Digital media watermarking is a popular image forensic problem. It requires to embed a signature into a sound, image, video,
3D mesh, etc.

An good source of information regarding digital watermarking is the book

One can also visit the BOWS-2 challenge homepage for a state of the art digital watermarking implementation.

We consider here a robust watermarking embedding, i.e. the goal is to embed a watermark that is both impercevable and difficult
to remove (by attack such as compression, denoising, adding noise, blurring, etc).

This is somehow conflicting goals since impercevable information is likely to be removed by an efficient compression or denoising
algorithm. An efficient watermarking scheme should thus use more clever tools than state of the art denoising/compression
algorithms.

Note also that we perform here "0 bit" watermarking, i.e. we do not embed a meaningful message within the watermarking. We
are only interested in testing the presence of a given watermark.

Here we bench a wavelet method for the embedding of a single watermark. We check how much the watermark can be detected after
various attack. Depending on a probability of false alarm, we compute the probability of detecting the watermark.

Watermark Embedding

A watermark is computed as a weighted random vector that is added to the wavelet coefficient.

The weighting of the watermark vector takes into account the amplitude of the host coefficient in order to reduce visual distortion.
This also increases the robustness to denoising and compression attacks.

The coefficients to be watermarked \(x_0 \in \RR^P \) is only a subset \( x_0 = (a_i)_{i \in I} \) of the total set of coefficients,
where \(\abs{I}=P\).

We select here only the fine scale wavelets.

A = ones(n); A(1:2^Jmin,1:2^Jmin) = 0;
I = find(A(:));
P = length(I);

Extract the coefficients \(x_0\).

x0 = a(I);

The watermarking is embedded using a multiplicative rule as \[ x_i = (x_0)_i + \rho \abs{ (x_0)_i } w_i \] where \(w\) is
a random Gaussian vector and where \(\rho > 0\) is a constant that ensure that \(\norm{x_0-x}\) is a given deviation value.

Generate the base watermark vector \(w \in \RR^P\).

w = randn(P,1);

Target embedding PSNR (should be quite large for the embedding to be unoticeable).

Watermark Detection

The watermark is detected (or not detected) from an input vector \(y \in \RR^P\) using a detector function \(C(y,w) \in \RR\)
where \(w \in \RR^P\) is the base watermark vector. Usually, a large value of \(C\) means that \(y\) is likely to come from
a watermarked content.

The detection is carried over by a simple thresholding, and the watermark is declared to be present if \[ C(y,w)>T \] where
\(T \in \RR\) is a threshold that should be set to guarantee a given probability of false alarms (i.e. ratio of contents
declared to be watermarked whereas they were not watermarked).

The detection corresponds to an hypothesis testing. One assumes that \(y=A(x)\) is obtained by attacking some vector \(x\),
and one has the following alternative depending on wether the content \(x\) is watermarked or not: \[ \choice{ (\Hh_0) \quad
x=x_0+\rho\abs{x_0}w, \\ (\Hh_1) \quad x=x_0. } \]

The two important quantities to monitor is the probability of false alarms \[ p_{\text{FA}} = \PP_w\pa{ C(y,w)>T \:\vert\:
\Hh_1 } \] and the probability of true positives \[ p_{\text{TP}} = \PP_w\pa{ C(y,w)>T \:\vert\: \Hh_0 }. \] Note that here
\(\PP_w\) refers to the probability of an event with respect to the randomization of \(w\).

The goal is to design a watermarking scheme (i.e. an embedding strategy and a detection strategy) in order to maximize \(p_{\text{TP}}\)
for a given \(p_{\text{FA}}\).

To estimate easily the probability of false alarm, we make the asumption that \(y\) is close enough to \(x_0\) to estimate
\(p_{\text{FA}}\) on the clean original signal \[ p_{\text{FA}} \approx \PP_w( C(x_0,w)>T ) \]

Exercice 3: (check the solution) Using a Monte Carlo simulation (generation of the order of \(10^3\) watermarks, display the histogram of the repartition
of \(C(x_0,w)\). Compute the variance \(\sigma_0^2\) of this distribution.

exo3;

We make another approximation : we approximate this density probability with a Gaussian density of mean 0 and variance \(\si_0^2\).
Under this assumption, one has \[ p_{\text{FA}} \approx 1 - G_{\si_0}(T) = 1 - \frac{1}{2} \pa{ 1 + \text{erf}\pa{\frac{T}{\sqrt{2}
\si_0}} } \] where \(G_{\si_0}\) is the cumulative density function of the Gaussian of variance \(\si_0^2\).

Hence one can use the threshold \[ T = \sqrt{2} \sigma_0 \text{erf}^{-1}(1-2 p_{\text{FA}}) \] This is an example of determination
of threshold \(T\) given a value of \( p_{\text{FA}} \).

pfa = 1e-3;
T = sqrt(2)/2 * sigma0 * erfinv(1-2*pfa);

Actually, it is possible to compute exactly this probability of false alarm as \[ p_{\text{FA}} = 1 - B(T^2 ; 1/2, (P-1)/2),
\] where \(B\) is the incomplete beta function (use betainc function) and \(P\) is the dimension.

Exercice 4: (check the solution) Compare, for various values of \(T\) the estimation obtained by the Gaussian approximation with the true value obtained
with the incomplete beta function.

exo4;

Quantization Attack

A compression attack is simulated by quantizing the wavelet coefficients. We consider here a dead zone quantization attack.

Quantization step \(\tau\) (the larger, the more aggressive the compression.

Exercice 5: (check the solution) Compute, by Monte Carlo sampling (i.e. draw at random many \(w\)) the distribution of \(C(A(x),w)\) for \(x = x_0 + \rho
\abs{x_0} w\). Store the different realization of \(C(A(x),w)\) in a vector c. Note: the value of \(\rho\) should be recomputed for each \(w\).

exo5;

Exercice 6: (check the solution) Compute, for a varying value of \( p_{\text{FA}} \), the corresponding value of \( p_{\text{TP}} \). Display the resulting
curve (ROC curve). This computation should be performed experimentally using e.g. 1000 random sampling.

exo6;

Exercice 7: (check the solution) Try different attack strengths, by changing the value of \(\tau\). For a \(p_{\text{FA}}=10^{-6}\), determine the value
of \(\tau\) for witch \(p_{\text{TP}}\) drops bellow \(0.2\).