function [SignalImgs,NoiseImgs,Nominal] = genFaces (N,SNR,SigPos,type,downscl)

%
% function [SignalImgs,NoiseImgs,Nominal] = genFaces(N,SNR,SigPos,downscl,type)
%
%  This function computes examples of cartoon face images together with noise
%  images. The positions of facial features are varied according to a Gaussian
%  density. White Gaussian pixel noise is added to all images according to the
%  desired signal-to-noise ratio (SNR).
%
%  Input Variables:
%
%       N               scalar          Number of images to generate. N noise and
%                                       N face mages will be generated.
%
%       SNR             scalar          Signal-to-noise ratio in dB. Try values 
%                                       between 1 and 20.
%
%       SigPos          scalar          Standard deviation of the Gaussian random 
%                                       variables that determine the positions of 
%                                       the facial features. 
%
%       donwscl         scalar          Indicates the number of times the images
%                                       should be donwsampled by a factor of 2.
%                                    Resulting image sizes for downscl = {0,1,2}
%                                       are {50x50, 25x25, 13x13}.
%
%       type           {1 or 2}         Indicates which face type to generate
%
%  Outputs Variables:
%
%       SignalImgs      NPix x N        Images containing faces stored as columns.
%
%       NoiseImgs       NPix x N        Images contianing just noisy pixels as columns.
%
%       Nominal         PixY x PixX     Nominal image with all features in the average
%                                       positions. This can be use for matched filtering 
%                                       of the entire face.
%
%  Image size is currently hard-wired to 50 x 50 pixels. If you make the positional
%  variance too large, the function might crash, since features move off the image.
%
%  (c) in 1999 by California Institute of Technology
%
%  Markus Weber
%


if (nargin < 3)
  type = 1;
end
  
% Call the script that defines how the parts look
faceParts;

NParts = length(Parts);

% Create two index vectors that will allow us to place the parts inside
% the large image matrix.

for p = 1 : NParts,

        % We assume that all sizes are odd
        partGridY{p} = [0 : size(Parts{p}, 1) - 1]' - floor(size(Parts{p}, 1) / 2);
        partGridX{p} = [0 : size(Parts{p}, 2) - 1]' - floor(size(Parts{p}, 2) / 2);
end

ImSize = [50 50]';
NPixels = prod(ImSize);

% mu and Sig correspond to points in order x1, x2, ..., xN, y1, y2, ..., yN
if (type == 1)
  mu = [ -4 4 0 0 -7 -7 0 9 ]';
else
  mu = [ -7 7 0 0 -5 -5 0 6 ]';
end        

 

Sig = eye(2 * NParts) * SigPos ^ 2;

mu = mu + [ones(NParts, 1) * ImSize(2) / 2; ones(NParts, 1) * ImSize(1) / 2];

% Create image where all parts are in nominal (mean) position

Nominal = zeros(ImSize(1), ImSize(2));

for p = 1 : NParts,

        Nominal(mu(p + NParts) + partGridY{p}, mu(p) + partGridX{p}) = ...
                Parts{p};
end

for i = 1 : downscl,

        Nominal = downsample(Nominal);
end

% Compute total energy in signal 
E = sum(sum(Nominal .^ 2));

% Use defition of SNR to compute pixels noise stddev needed for desired SNR
SigNoise = sqrt(E) / (10 ^ (SNR / 20));

% Draw random variables for part positions in all images
Pos = round(randvec(N, mu, Sig));

SignalImgs = zeros(length(Nominal(:)), N);

for i = 1 : N,

        img = zeros(ImSize(1), ImSize(2));

        for p = 1 : NParts,

        img(Pos(p + NParts, i) + partGridY{p}, Pos(p, i) + partGridX{p}) = ...
                        Parts{p};

        end

        img = img + SigNoise * randn(size(img));

        %imagen(img); drawnow;

        for j = 1 : downscl,
                
                img = downsample(img);
        end

        SignalImgs(:, i) = img(:);
end

NoiseImgs = zeros(size(SignalImgs, 1), N);

for i = 1 : N,

        img = SigNoise * randn(ImSize(1), ImSize(2));

        for j = 1 : downscl,
                img = downsample(img);
        end
        
        NoiseImgs(:, i) = img(:);
end