function [imw,code] = watermark(im,signature,strength,density);

% [imw,code] = WATERMARK(inputfile,signature,strength,density);
%
% Watermarks the image IM and returns the watermarked
% image IMW and the retrieval code.
%
% signature:  Watermark, STRING
% strength:   Signature Strength, 0..1
% density:    Signature Density, 0..1
%
% e.g.: [imw,code]=watermark(im,'signature',0.1,0.55);
%
%
% University of Rhode Island
% Dept. of Electrical Engineering
% Roger Meier, May '99


RGB=double(im);

% Calculating luminance of the image file
% ---------------------------------------
lum=0.299*RGB(:,:,1) + 0.587*RGB(:,:,2) + 0.114*RGB(:,:,3);
		
% Calculating seed
% ----------------
rand('state',sum(100*clock));
seed=round(9999*rand);
      
% Feeding random number generator with seed
% -----------------------------------------
rand('state',seed);

% Calculating embedding positions
% -------------------------------
H=size(RGB,1);
W=size(RGB,2);
emb = (rand(H,W) < density);
[y,x] = find(emb);
pos = [y x];
posnum = length(x);

% Converting signature into bitstring
% -----------------------------------
bit1=dec2bin(double(signature),7);
bit2='';
for k = 1:size(bit1,1)
	bit2=[bit2 bit1(k,:)];
end

% Adding the bits '0' and '1' to the bitstring
% to allow to define an adaptive threshold which
% improves signature retrieval. This also defines
% a geometrical reference which is used to counter
% geometrical attacks.
% ------------------------------------------------
bit2=['0' '1' bit2];

% Converting the bitstring to a real vector
% -----------------------------------------
bitlength=length(bit2);
bit=zeros(1,bitlength);
for k = 1:bitlength
	bit(k)=str2num(bit2(k));
end

% Feeding random number generator with seed
% -----------------------------------------
rand('state',seed);

% 'Shuffeling' the positions of the embedding bits
% ------------------------------------------------
order=randperm(posnum);
pos=pos(order,:);

% Generating the watermark bitstring of length fig.posnum
% by lining up copies of the bitstring
% -------------------------------------------------------
Watermark=zeros(posnum,1);
for k=1:floor(posnum/bitlength)
   Watermark(bitlength*(k-1)+1:bitlength*k)=bit';
end
remain=posnum-bitlength*floor(posnum/bitlength);
Watermark(posnum-remain+1:posnum)=bit(1:remain)';

% Watermarking the image
% ----------------------
B=RGB(:,:,3);
Bnew=B;
y=pos(:,1);
x=pos(:,2);
num=y+H*(x-1); % Converting XY Position to Pixel number
Bnew(num) = B(num) + ((2*Watermark-1) .* lum(num) .* strength);

% Creating new image
% ------------------
under=Bnew < 0;  % setting all elements < 0 in new blue channel to zero
under=find(under);
Bnew(under)=0;
over=Bnew > 255; % setting all elements > 255 in new blue channel to 255
over=find(over);
Bnew(over)=255;
RGBnew=RGB;
RGBnew(:,:,3)=Bnew;

imw=uint8(RGBnew);

% Calculating retrieval code
% --------------------------
if seed > 1000
   seedcode=num2str(seed);
elseif seed > 100
   seedcode=['0' num2str(seed)];
elseif seed > 10
   seedcode=['00' num2str(seed)];
else
   seedcode=['000' num2str(seed)];
end
if density == 1
   denscode=num2str(100*density);
else
   if density >= 0.1
      denscode=['0' num2str(100*density)];
   else
      denscode=['00' num2str(100*density)];
   end
end
lengthcode=num2str(bitlength);

code=[seedcode denscode lengthcode];