%
% Matlab function to demonstrate
% the spread-spectrum technique for
% watermarking. Refer project webpage
% section on this for the details
%
% To be consistent across insertion and extraction
% process I decided on the following parameters
%
% Return values
% -------------
% D   - the output image matrix. 
% V   - the original DCT coefficients
% X   - ~N(0,1) random vector of length n ; used in similarity function
% I   - the indices of the high value DCT coefficients
% Input parameters:
% -----------------
% Din - the image matrix. (usually got from the imread of
%     an image GIF)
% n - the number of watermark bits to embed (say, 1000)
% mode - the demo mode you want to test. possible values are
%   'random' which is a test for the watermark detector
%            of original watermark vs random watermarks. the
%            output D matrix is the same as Din matrix
%   'jpeg'   which is used as a test for JPEG compression
%            and generates a return image matrix after writing/
%            reading to a JPEG file
%   'scaledjpeg' which is used to test the effect of both 
%                scaling and jpeg compression
%    
function [D, V, X, I] = spread_watermark_insert(Din, n, mode)

% convert to double
mat = double(Din);

% perform DCT on the whole image
% using blockproc
fun = @dct2;
dct=blkproc(mat,[8 8], fun);

% find the top n dct coefficients 
mask = findnmax(dct, n);

% get the correct indices for the top
% n coefficients
[a,b] = find(mask==1);

% for the number of watermark bits
for i=1:n,
    
    % genarate a random number from 
    % the randn function (mean = 0, std = 1)    
    X(1,i) = randn(1);
    
    % store the original values of the DCT
    V(1,i) = dct(a(i), b(i));
    
    % replace the original values with the 
    % following invertible formula
    % v(i) = v(i)* (1 + alpha* x(i))
    % here alpha is the scaling coefficent and
    % as in the paper 0.1 is chosen
    alpha = 0.1;    
    dct(a(i), b(i)) = dct(a(i), b(i)) * (1 + (alpha*X(1,i)) );     
    
end

% go back to the spatial domain
idct=blkproc(dct, [8 8], 'idct2'); 
mat=uint8(idct); 


% get the indices of the high DCT coefficients
I = [a b];

% generate the appropriate image matrix
% D based on the mode of the operation
% (quality numbers assumed randomly. can easily
% add another parameter for the function for quality
mode
if strcmp(mode, 'jpeg'),
    display('writing out to jpeg');
    imwrite(mat,'tmpstego.jpg','JPEG', 'Quality', 10, 'Mode', 'lossy');
    D = double(imread('tmpstego.jpg', 'JPEG'));
else,
    if strcmp(mode, 'scaledjpeg'),
        display('writing out to scaled jpg');
        imwrite(mat, 'tmpstego.jpg', 'JPEG', 'Quality',90);
        D = double(imread('tmpstego.jpg', 'JPEG'));
    else,
        D = double(Din);
    end
end        