Image processing Sum Square Diffrence

I'm having trouble getting my code to work in Matlab. I try to write code witch use the below formula

I have to find in image (I) template (T) with other dimensions

My actual code

clear;clc;

I=imread('00.bmp');
I=im2double(I);
[W H]=size(I);

T=imread('wzorp.png');
T=rgb2gray(T);
T=im2double(T);
[w h]=size(T);


for x=1:W-w+1
    for y=1:H-h+1
        for x1=1:w-1
            for y1=1:h-1
                R(x,y)=(T(x1,y1)-I(x+x1,y+y1))^2
            end
        end
    end
end

Answer 1

[Thiago Henrique Gomes Lobato]

What exactly is the problem you're having? If the idea is just o find a template in the image you can probably find better results using image correlation, as example the normxcorr2 function.

[Oskar Mencwel]

the problem is the very long running time of the current code, i want to compare the results of use cross correlation, correlation coefficient and suared dufference, but I'm not sure if I wrote the function of suared dufference code correctly

[Thiago Henrique Gomes Lobato]

You have many loops, which are very unefficient in Matlab, that's why it is so slow. One can partially solve the problem adding an extra vectorization as follows (I changed the image):

clear;clc;
I=imread('cameraman.tif');
I=im2double(I);
I = padarray(I,[52,52],0);

[W H]=size(I);
T=imread('cameraman.tif');
%T=rgb2gray(T);
T=im2double(T);
[w h]=size(T);
R = zeros(W-w+1,H-h+1);
for x=1:W-w+1
    for y=1:H-h+1
        R(x,y) = sum( sum( (( T-I(x:x+w-1,y:y+h-1) ).^2)  ));
    end
end

figure,imshow(1-R/max(R(:))) % 1 means best template match, which coincides here in the center of R

You can write the loop as a mex-file to find a speed-up. The equation is right, the only limitation in it is that you don't do any padding in your image, so template in corners would probably not be founded, you can solve this by just padding your original image with the padarray I wrote in the example. If you also want to have the matrix R of the same size as the image I, just make the pad to compensate your template size.

[Image Analyst]

This doesn't look like a correlation, which is a multiplication and then a sum.

[Thiago Henrique Gomes Lobato]

Is it indeed just a difference metric, a proper correlation should work in almost every case better.

Answer 2

[Image Analyst]

Try this code based on convolution. It's fast = 0.16 seconds on my computer:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures if you have the Image Processing Toolbox.
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 22;
% Check that user has the Image Processing Toolbox installed.
hasIPT = license('test', 'image_toolbox');
if ~hasIPT
  % User does not have the toolbox installed.
  message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
  reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
  if strcmpi(reply, 'No')
    % User said No, so exit.
    return;
  end
end
%===============================================================================
% Read in a demo image.
folder = fileparts(which('peppers.png')); % Determine where demo folder is (works with all versions).
baseFileName = 'eight.tif';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
  % Didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
grayImage1 = imread(fullFileName);
% Get the dimensions of the image.  numberOfColorBands should be = 3.
[rows, columns, numberOfColorBands] = size(grayImage1);
if numberOfColorBands > 1
	fprintf('This image is RGB.  I will change it to gray scale.\n');
	grayImage1 = grayImage1(:, :, 2);
end
% Display the original image.
subplot(2, 2, 1);
imshow(grayImage1);
axis('on', 'image');
title('Image 1', 'FontSize', fontSize);
impixelinfo;
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'Outerposition', [0, 0.1, 1, 0.9]);

% Now get gray image #2
baseFileName = 'coins.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
  % Didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
grayImage2 = imread(fullFileName);
% Get the dimensions of the image.  numberOfColorBands should be = 3.
[rows, columns, numberOfColorBands] = size(grayImage2);
if numberOfColorBands > 1
	fprintf('This image is RGB.  I will change it to gray scale.\n');
	grayImage2 = grayImage2(:, :, 2);
end
% Make the images the same size
grayImage2 = imresize(grayImage2, size(grayImage1));

% Display the original image.
subplot(2, 2, 2);
imshow(grayImage2);
axis('on', 'image');
title('Image 2', 'FontSize', fontSize);
impixelinfo;

tic
% Create a kernel
windowWidth = 3;
kernel = ones(windowWidth, windowWidth)
kernel(2, 2) = 0;
% Convolve the second image with the kernel.
outputImage = conv2(double(grayImage2), kernel, 'same');
% Now 
subplot(2, 2, 3);
imshow(outputImage, []);
axis('on', 'image');
title('Convolution Image based on image 2', 'FontSize', fontSize);

% Now subtract and square.
scalingFactor = windowWidth^2 - 1;
imageR = (scalingFactor * double(grayImage1) - outputImage) .^ 2;
% Display the original image.
subplot(2, 2, 4);
imshow(imageR, []);
axis('on', 'image');
title('R Image', 'FontSize', fontSize);
impixelinfo;
toc