function [x, y]=SSDA_Match(img, subimg)
%-------------------------------------------------------------------------%
%该函数用来查找原始图像img中的子区块subimg的位置，位置值取匹配区域的左上角坐标。
%函数使用互相关法实现。为了减小搜索区域，先将原图像和子图像抽取至较小范围，获取
%相应位置后，再在原图像相应点附近详细搜索。
%返回值x、y是相应的坐标值。
%-------------------------------------------------------------------------%
[M, N] = size(img); %获取原图像、子图像大小
[Msub, Nsub] = size(subimg);
Mtmp = Msub; Ntmp = Nsub;
img = double(img);
subimg = double(subimg);
if (M > Msub) & (N > Nsub)
times = 0;
while (Mtmp > 32) | (Ntmp > 32) %获取可压缩的次数times
Mtmp = Mtmp/2 ;
Ntmp = Ntmp/2 ;
times = times + 1;
end
img1 = img(1:2.^times:M, 1:2.^times:N); %得到压缩图像
subimg1 = subimg(1:2.^times:Msub, 1:2.^times:Nsub);
[Ms, Ns] = size(img1);
[Mssub, Nssub] = size(subimg1);
CorRes = zeros(Ms-Mssub+1, Ns-Nssub+1); %创建并且计算压缩图互相关结果
for delta_i = 1:Ms-Mssub+1
for delta_j = 1:Ns-Nssub+1
tmp1 = sum(sum(img1(delta_i:delta_i+Mssub-1, delta_j:delta_j+Nssub-1).*subimg1(1:Mssub, 1:Nssub)));
tmp2 = sqrt(sum(sum(img1(delta_i:delta_i+Mssub-1, delta_j:delta_j+Nssub-1).^2)));
tmp3 = sqrt(sum(sum(subimg1(1:Mssub, 1:Nssub).^2)));
CorRes(delta_i, delta_j) = tmp1/(tmp2*tmp3);
end
end
[tmp, k] = max(CorRes);
[tmp, kk] = max(tmp);
xs = k(kk);
ys = kk; %取得最大互相关坐标
x1 = 1 + (xs-1) .* (2.^times); %找到原始图像带搜索的区域
x2 = min([xs .* (2.^times) M]); %防止坐标越界
y1 = 1 + (ys-1) .* (2.^times);
y2 = min([ys .* (2.^times) N]);
CorRes1 = zeros(x2-x1+1, y2-y1+1); %创建并且计算原图互相关结果
for delta_i = x1:x2
for delta_j = y1:y2
tmp1 = sum(sum(img(delta_i:delta_i+Msub-1, delta_j:delta_j+Nsub-1).*subimg(1:Msub, 1:Nsub)));
tmp2 = sqrt(sum(sum(img(delta_i:delta_i+Msub-1, delta_j:delta_j+Nsub-1).^2)));
tmp3 = sqrt(sum(sum(subimg(1:Msub, 1:Nsub).^2)));
CorRes1(delta_i-x1+1, delta_j-y1+1) = tmp1/(tmp2*tmp3);
end
end
[tmp, k] = max(CorRes1);
[tmp, kk] = max(tmp);
delta_x = k(kk) - 1;
delta_y = kk - 1; %取得最大互相关坐标
x = x1 + delta_x; %得到精确坐标
y = y1 + delta_y;
% str1 = strcat('计算得到的精确匹配坐标为 [',num2str(x),','num2str(y),']')
else
error('子图的维数比待查图像大'); %子图大小越界
end
