SSD演算法

        誤差平方和演算法（Sum of Squared Differences，簡稱SSD演算法），也叫差方和演算法。實際上，SSD演算法與SAD演算法如出一轍，只是其相似度測量公式有一點改動（計算的是子圖與模板圖的L2距離）。這裡不再贅述。

matlab 程式碼：

<pre name="code" class="cpp">%%
%絕對誤差和演算法（SAD）
clear all;
close all;
%%
src=imread('lena.jpg');
[a b d]=size(src);
if d==3
    src=rgb2gray(src);
end
mask=imread('lena_mask.jpg');
[m n d]=size(mask);
if d==3
    mask=rgb2gray(mask);
end
%%
N=n;%模板尺寸，預設模板為正方形
M=a;%代搜尋影象尺寸，預設搜尋影象為正方形
%%
dst=zeros(M-N,M-N);
for i=1:M-N         %行
    for j=1:M-N
        temp=src(i:i+N-1,j:j+N-1);
        dst(i,j)=dst(i,j)+sum(sum((temp-mask).^2));
    end
end
abs_min=min(min(dst));
[x,y]=find(dst==abs_min);
figure;
imshow(mask);title('模板');
figure;
imshow(src);
hold on;
rectangle('position',[x,y,N-1,N-1],'edgecolor','r');
hold off;title('搜尋圖');
 

或者：

利用頻域的變化

function [I_SSD,I_NCC,Idata]=template_matching(T,I,IdataIn)
% TEMPLATE_MATCHING is a cpu efficient function which calculates matching 
% score images between template and an (color) 2D or 3D image.
% It calculates:
% - The sum of squared difference (SSD Block Matching), robust template
%   matching.
% - The normalized cross correlation (NCC), independent of illumination,
%   only dependent on texture
% The user can combine the two images, to get template matching which
% works robust with his application. 
% Both measures are implemented using FFT based correlation.
%
%   [I_SSD,I_NCC,Idata]=template_matching(T,I,Idata)
%
% inputs,
%   T : Image Template, can be grayscale or color 2D or 3D.
%   I : Color image, can be grayscale or color 2D or 3D.
%  (optional)
%   Idata : Storage of temporary variables from the image I, to allow 
%           faster search for multiple templates in the same image.
%
% outputs,
%   I_SSD: The sum of squared difference 2D/3D image. The SSD sign is
%          reversed and normalized to range [0 1] 
%   I_NCC: The normalized cross correlation 2D/3D image. The values
%          range between 0 and 1
%   Idata : Storage of temporary variables from the image I, to allow 
%           faster search for multiple templates in the same image.
%
% Example 2D,
%   % Find maximum response
%    I = im2double(imread('lena.jpg'));
%   % Template of Eye Lena
%    T=I(124:140,124:140,:);
%   % Calculate SSD and NCC between Template and Image
%    [I_SSD,I_NCC]=template_matching(T,I);
%   % Find maximum correspondence in I_SDD image
%    [x,y]=find(I_SSD==max(I_SSD(:)));
%   % Show result
%    figure, 
%    subplot(2,2,1), imshow(I); hold on; plot(y,x,'r*'); title('Result')
%    subplot(2,2,2), imshow(T); title('The eye template');
%    subplot(2,2,3), imshow(I_SSD); title('SSD Matching');
%    subplot(2,2,4), imshow(I_NCC); title('Normalized-CC');
%
% Example 3D,
%   % Make some random data
%    I=rand(50,60,50);
%   % Get a small volume as template
%    T=I(20:30,20:30,20:30);
%   % Calculate SDD between template and image
%    I_SSD=template_matching(T,I);
%   % Find maximum correspondence
%    [x,y,z]=ind2sub(size(I_SSD),find(I_SSD==max(I_SSD(:))));
%    disp(x);
%    disp(y);
%    disp(z);
%
% Function is written by D.Kroon University of Twente (February 2011)

if(nargin<3), IdataIn=[]; end

% Convert images to double
T=double(T); I=double(I);
if(size(T,3)==3) 
    % Color Image detected
    [I_SSD,I_NCC,Idata]=template_matching_color(T,I,IdataIn);
else
    % Grayscale image or 3D volume
    [I_SSD,I_NCC,Idata]=template_matching_gray(T,I,IdataIn);
end

function [I_SSD,I_NCC,Idata]=template_matching_color(T,I,IdataIn)
if(isempty(IdataIn)), IdataIn.r=[];  IdataIn.g=[]; IdataIn.b=[];  end
% Splite color image, and do template matching on R,G and B image
[I_SSD_R,I_NCC_R,Idata.r]=template_matching_gray(T(:,:,1),I(:,:,1),IdataIn.r);
[I_SSD_G,I_NCC_G,Idata.g]=template_matching_gray(T(:,:,2),I(:,:,2),IdataIn.g);
[I_SSD_B,I_NCC_B,Idata.b]=template_matching_gray(T(:,:,3),I(:,:,3),IdataIn.b);
% Combine the results
I_SSD=(I_SSD_R+I_SSD_G+I_SSD_B)/3;
I_NCC=(I_NCC_R+I_NCC_G+I_NCC_B)/3;

    
function [I_SSD,I_NCC,Idata]=template_matching_gray(T,I,IdataIn)
% Calculate correlation output size  = input size + padding template
T_size = size(T); I_size = size(I);%模板和影象的大小
outsize = I_size + T_size-1;%這個是擴充套件邊界後的大小

% calculate correlation in frequency domain
if(length(T_size)==2)
    FT = fft2(rot90(T,2),outsize(1),outsize(2));%計算模板的的頻域
    if(isempty(IdataIn))
        Idata.FI = fft2(I,outsize(1),outsize(2));%計算影象的頻域
    else
        Idata.FI=IdataIn.FI;
    end
    Icorr = real(ifft2(Idata.FI.* FT));%進行反變換 時域卷積，頻域相乘   得到二者的相關性
else
    FT = fftn(rot90_3D(T),outsize);
    FI = fftn(I,outsize);
    Icorr = real(ifftn(FI.* FT));
end

% Calculate Local Quadratic sum of Image and Template
if(isempty(IdataIn))
    Idata.LocalQSumI= local_sum(I.*I,T_size);%計算影象的區域性的差的平方
else
    Idata.LocalQSumI=IdataIn.LocalQSumI;
end

QSumT = sum(T(:).^2);%模板中全部畫素的和平方

% SSD between template and image
I_SSD=Idata.LocalQSumI+QSumT-2*Icorr;

% Normalize to range 0..1
I_SSD=I_SSD-min(I_SSD(:)); 
I_SSD=1-(I_SSD./max(I_SSD(:)));

% Remove padding
I_SSD=unpadarray(I_SSD,size(I));

if (nargout>1)
    % Normalized cross correlation STD
    if(isempty(IdataIn))
        Idata.LocalSumI= local_sum(I,T_size);%計算出影象中模板快內的畫素之和
    else
        Idata.LocalSumI=IdataIn.LocalSumI;
    end
    
    % Standard deviation
    if(isempty(IdataIn))
        Idata.stdI=sqrt(max(Idata.LocalQSumI-(Idata.LocalSumI.^2)/numel(T),0) );%numl元素的個數
    else
        Idata.stdI=IdataIn.stdI;
    end
    stdT=sqrt(numel(T)-1)*std(T(:));
    % Mean compensation
    meanIT=Idata.LocalSumI*sum(T(:))/numel(T);
    I_NCC= 0.5+(Icorr-meanIT)./ (2*stdT*max(Idata.stdI,stdT/1e5));

    % Remove padding
    I_NCC=unpadarray(I_NCC,size(I));
end

function T=rot90_3D(T)
T=flipdim(flipdim(flipdim(T,1),2),3);

function B=unpadarray(A,Bsize)
Bstart=ceil((size(A)-Bsize)/2)+1;
Bend=Bstart+Bsize-1;
if(ndims(A)==2)
    B=A(Bstart(1):Bend(1),Bstart(2):Bend(2));
elseif(ndims(A)==3)
    B=A(Bstart(1):Bend(1),Bstart(2):Bend(2),Bstart(3):Bend(3));
end
    
function local_sum_I= local_sum(I,T_size)
% Add padding to the image
B = padarray(I,T_size);%對影象進行填充

% Calculate for each pixel the sum of the region around it,
% with the region the size of the template.  計算模板區域內的和  就是SSD
if(length(T_size)==2)
    % 2D localsum
    s = cumsum(B,1);
    c = s(1+T_size(1):end-1,:)-s(1:end-T_size(1)-1,:);
    s = cumsum(c,2);
    local_sum_I= s(:,1+T_size(2):end-1)-s(:,1:end-T_size(2)-1);
else
    % 3D Localsum
    s = cumsum(B,1);
    c = s(1+T_size(1):end-1,:,:)-s(1:end-T_size(1)-1,:,:);
    s = cumsum(c,2);
    c = s(:,1+T_size(2):end-1,:)-s(:,1:end-T_size(2)-1,:);
    s = cumsum(c,3);
    local_sum_I  = s(:,:,1+T_size(3):end-1)-s(:,:,1:end-T_size(3)-1);
end