function varargout=retrieve(imw,code,crosssize);

% signature = RETRIEVE(im,code,crosssize);
%
% 		Retrieves the watermark embedded in IMW
% 		and returns the signature.
%
% 		code      : Retrieval code, STRING
% 		crosssize : Size of crossize neighborhood, INTEGER
%
% 		e.g.: signature = retrieve(imw,'241805579',3);
%
%
% [signature,delta] = retrieve(imw,code,crosssize);
% 		also returns the averaged deltas
%
%
% University of Rhode Island
% Dept. of Electrical Engineering
% Roger Meier, May '99

RGB=double(imw);

% Decoding code
% -------------
seed=str2num(code(1:4));
density=str2num(code(5:7))/100;
bitlength=str2num(code(8:end));

% 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);

% Extracting blue channel
% -----------------------
B=RGB(:,:,3);

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

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

% Estimating original blue values
% using cross shaped neighborhood
% -------------------------------
cross=zeros(2*crosssize+1,2*crosssize+1);
cross(:,crosssize+1)=1;
cross(crosssize+1,:)=1;
cross(crosssize+1,crosssize+1)=0;
Bold=1/(4*crosssize) * filter2(cross,B,'same');

% Calculating the differences
% ---------------------------
delta=B-Bold;

% Extracting deltas at embedding positions
% ----------------------------------------
y=pos(:,1);
x=pos(:,2);
num=y+H*(x-1); % Converting XY Position to Pixel number
delta=delta(num);

% Extracting and averaging signature bits
% ---------------------------------------
bit=zeros(bitlength,1);
for k=1:floor(posnum/bitlength)
   bit=bit+delta(bitlength*(k-1)+1:bitlength*k);
end
bit=bit/floor(posnum/bitlength);

% Calculate adaptive treshhold
% ----------------------------
thres=(bit(1)+bit(2))/2;

% Retrieving signature
% --------------------
valid = (bitlength-2)/7 == round((bitlength-2)/7);
if valid
	sigbin = bit > thres;
	sigbin = sigbin(3:end);
	sigasc = zeros(1,round(length(sigbin)/7));
	for k = 1:length(sigasc)
   	bin = '';
   	for l = 1:7
      	bin=[bin num2str(sigbin(l+7*(k-1)))];
   	end
   	sigasc(k)=bin2dec(bin);
	end
   signature=char(sigasc);
else
   signature='INVALID RETRIEVAL CODE !!!';
end


switch nargout
case 0,
   varargout = {signature};
case 1,
   varargout = {signature};
case 2,
   varargout = {signature,bit};
end
