function [J, W] = CoxWmk(I, alpha, N, W)
%COXWMK Insert a random watermark.
%   [J, W] = COXWMK(I) watermarks the grayscale picture I using
%   Cox's scheme [2] with default parameters. J is the watermarked
%   picture and W is the watermark.
%   The watermark W is drawn from a Normal random process N(0,1)
%   if it is not given by the user as last parameter.
% 
%   [J, W] = COXWMK(I, Alpha, N) watermarks the grayscale picture I
%   using modifying the Nth largest DCT coefficients v_i in the
%   following way: v_i := v_i * (1 + Alpha * w_i) where w_i
%   are the components of the watermak.
%
%   Use IMSHOW(mat2gray(J)) to display the output.
%
%   See: COXDETECT

%   Fabien A.P. Petitcolas 26-11-97
%   Copyright (c) 1997 by the University of Cambridge
%   $Revision: 0.6$

%   References:
%        1) Ross J. Anderson, editor. Information hiding: first
%           international workshop, volume 1174 of Lecture notes
%           in computer science. University of Cambridge, Isaac
%           Newton Institute, Springer Verlag, Berlin, Germany,
%           May 1996.
%        2) Ingemar J. Cox, Joe Kilian, Tom Leighton, and Talal
%           Shamoon. A secure, robust watermark for multimedia.
%           In Anderson [1], pages 183-206

if (nargin == 1)
   N = 1000;
   alpha = 0.1;
   W = randn(1,N);
end
if (nargin == 2)
   N = 1000;
   W = randn(1,N);
end
if (nargin == 3)
   W = randn(1,N);
end

sI = size(I);
if ((sI(1) * sI(2)) < N)
   error('Image too small or too many coefficients.');
end

% Generate the random watermark
%W = randn(1,N);

% Compute the DCT of the image
DCTI = dct2(I);

% Find the N largest coefficients in the DCT matrix
% Better if the extraction of the N largest was done
% at the same time than the computation of the DCT...
Index = FindNLargest(abs(DCTI), N);

% Modify these coefficients
for i = 1:N
   DCTI(Index(1,i),Index(2,i)) = DCTI(Index(1,i),Index(2,i)) * (1 + alpha * W(i));
end

% Simply take the inverse DCT of the modifyied matrix
J = idct2(DCTI);