clear all;

% read in all the images and square them
image1 = double(imread('f5.png'));
image1sq = image1.^2;

imageNameList = {'f1.png', 'f2.png', 'f3.png', 'f4.png', 'f6.png', 'f7.png', 'f8.png', 'f9.png'};
imageStack = zeros(size(image1,1), size(image1,2), size(imageNameList, 2));

for i=1:size(imageNameList, 2)
image = double(imread(imageNameList{i}));
imageStack(:,:,i) = image;
end

figure('Name', 'Final Image', 'Position', [0 0 512 512]);
finalImage = zeros(1024, 1024);
finalImage(257:512, 257:512) = image1;
imagesc(finalImage); axis square; colormap gray; hold on


for h=1:size(imageStack, 3)

    image2 = imageStack(:,:,h);
    image2sq = image2.^2;
    
for i=1:size(image1,1) - 32:256
    for j=1:size(image1,2) - 32:256
    % select a sub image from the first image read
            subImage1 = image1(i:(i + 31), j:(j+31));

            % find the sum total squared of the sub image for matching
            totsq = sum(sum(subImage1.^2));

            %pad all the images
            s1_1 = size(image1,1);
            s2_1 = size(image1,2);

            subImagePad = zeros(s1_1, s2_1);
            subImagePad(1:32, 1:32) = subImage1;

            mask = zeros(s1_1, s2_1);
            mask(1:32, 1:32) = 1.0;

            image1F = fft2(image1);
            image2F = fft2(image2);

            subImgF = fft2(subImagePad);

            image1sqF = fft2(image1sq);
            image2sqF = fft2(image2sq);

            maskF = fft2(mask);

            cc1F = image1F.*conj(subImgF);
            cc2F = image2F.*conj(subImgF);

            cc1 = abs(ifft2(cc1F));
            cc2 = abs(ifft2(cc2F));

            cct1 = abs(ifft2(image1sqF.*conj(maskF))) - 2 * cc1 + totsq;
            cct2 = abs(ifft2(image2sqF.*conj(maskF))) - 2 * cc2 + totsq;

            [m1, indx1] = min(cct1(:))
            [m2, indx2] = min(cct2(:))
            if m2 < 1.0e2
                disp('*******************************Found Match');
                break
            end
    end %for j
    if m2 < 1.0e2
        break
    end
end %for i

    if m2 > 1.0e2
        continue
    end %for h


[xoff1, yoff1] = ind2sub(size(cct1), indx1)
[xoff2, yoff2] = ind2sub(size(cct2), indx2)


xoff1 = xoff1 + 256;
yoff1 = yoff1 + 256;
finalImage(xoff1 - xoff2:255 + xoff1 - xoff2, yoff1 - yoff2:255 + yoff1 - yoff2) = image2;
figure(1);
imagesc(finalImage); axis square; colormap gray; hold on

end


%imagesc(finalImage); axis square; colormap gray; hold on