Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This numerical tour explores robust watermarking over the wavelet domain.
Many thanks to Patrick Bas and Teddy Furon for their useful advices on digital image watermarking.
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('solutions/coding_5_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
Digital Watermarking and Steganography, 2nd Ed, Ingemar J. Cox, Matthew L. Miller, Jeffrey Bloom, Morgan Kaufmann, 2007.
One can also visit the <http://bows2.gipsa-lab.inpg.fr/ 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.
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.
Load an image $f \in \RR^N$ of $N = n \times n$ pixels.
n = 256;
name = 'hibiscus';
f = load_image(name, n);
f = rescale(sum(f,3));
Display the original image.
clf;
imageplot(f);
Shortcut for the wavelet transform $\Psi : f \in \RR^N \mapsto a \in \RR^N$ that maps an image $f$ to wavelet coefficients $a$. We note its inverse $\Psi^{-1}$ using the shortcut |PsiS|.
Jmin = log2(n)-2;
Psi = @(f)perform_wavelet_transf(f, Jmin, +1);
PsiS = @(a)perform_wavelet_transf(a, Jmin, -1);
Compute the wavelet coefficients.
a = Psi(f);
Display the wavelet coefficients.
clf;
plot_wavelet(a,Jmin);
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).
psnr_embedding = 50;
Exercise 1
Compute |rho| so that |PSNR(y,x0,1)=snr_embedding|.
exo1()
rho = 0.0795.
%% Insert your code here.
Exercise 2
According to you, for which PSNR the watermark becomes unoticeable?
exo2()
%% Insert your code here.
Perform the embedding $x=x_0+\rho\abs{x_0}w$.
x = x0 + rho*abs(x0).*w;
The distortion of the embedding is measured using the PSNR $$ \text{PSNR}(x,x0) = -20 \log_{10}( \norm{x-x0}/\sqrt{P} ). $$
Check the PSNR of embedding.
disp(['PSNR(x,x0) = ' num2str(psnr(x,x0,1), 3) 'dB.']);
PSNR(x,x0) = 50dB.
Given the watermarked coefficients $x \in \RR^P$, a watermarked image $f_1 \in \RR^N$ is reconstructed using the inverse wavelet transform $\Psi^{-1}$ as $$ f_1 = \Psi^{-1}(a_1) \qwhereq (a_1)_i = \choice{ x_i \qifq i \in I, \\ a_i \quad\text{otherwise.} }$$
Compute the image with the watermark embedded.
a1 = a; a1(I) = x;
f1 = PsiS(a1);
Display the watermark $ \delta = \Psi^{-1}(a-a_1) = f - f_1 $ over the spacial domain (with contrast boosting).
delta = f-f1;
clf;
imageplot( clamp(delta/std(delta(:)),-3,3) );