function I = GrowImage(S, n, w, h)
if n/2 == round(n/2)
    error('n must be odd');
end
nh = (n-1)/2;
Ipad = init(S,n,w,h,nh);
S = double(S);
PixelList = GetUnfilledNeighbors(Ipad);
while ~isempty(PixelList)
    for i = 1:size(PixelList,2)
        Pixel = PixelList(:,i);
        Template = GetWindow(Pixel, Ipad, nh);
        % Template may contain -1 or -2; we don't use those pixels.
        BMs = FindMatches(Template, S);
        BM = RandomPick(BMs);
        Ipad(Pixel(1),Pixel(2)) = BM;
    end
    PixelList = GetUnfilledNeighbors(Ipad);
end
I = Ipad((nh+1):(h+nh), (nh+1):(w+nh));  

function BM = RandomPick(BMs)
% Pick a column at random


function BMs = FindMatches(Tem, S)
% This function will call SSD to compute the SSD between the template and every
% equal-sized window in S.  If the Template contains a -1, we want to
% ignore that pixel in computing the SSD.  Then we want to return the
% pixel values in S at the center of every region that has a distance 
% to the template that is no more than 1.1 times the distance of the 
% nearest region.  So BMs will be a vector of pixel intensities.

  
function D = SSD(S, Tem, M)
% D will be the same size as S.  It will give the SSD between Tem and the
% neighborhood around each pixel in S.  M is a mask, the size of Tem, that
% contains 0s and 1s.  1 means a pixel in the Template is valid and should
% be used in computing SSD, 0 means don't use that pixel (it's -1 or -2).

function Tem = GetWindow(P, I, nh);
% Extract an nhxnh window from I that is centered on pixel P.

function PL = GetUnfilledNeighbors(I)
% We want to find every -1 that is next to something that's not -1.
% For speed, it is important to do this without looping, using imfilter and
% find.   
% These are returned as a 2xn matrix, where every column gives the position
% of a pixel that is -1, but next to something not -1 (or -2).  The first row of PL
% gives the row position of the pixels, the second row gives the column
% position.

function Ipad = init(S,n,w,h,nh)
% We will pad the image with -2s so that we don't have to check boundary
% conditions when we index into it.  
% -1 will be a pixel we still need to fill in.
% We will start out Ipad with a 3x3 section in the center that we take
% randomly from the sample, S.
% Note: this function is complete as is; you don't need to modify it.
Ipad = [-2*ones(h+2*nh,nh), [-2*ones(nh,w+nh); ...
            [-1*ones(h,w), -2*ones(h,nh)]; -2*ones(nh,w+nh)]];
[Sh,Sw] = size(S);
r = 1+floor((Sh-2)*rand);
c = 1+floor((Sw-2)*rand);
hc = round(nh+h/2);
wc = round(nh+w/2);
Ipad(hc:(hc+2), wc:(wc+2)) = S(r:(r+2), c:(c+2));