clc;
clear all;
close all;
A = imread('bd1.JPG');%读取无火焰图像
B = imread('bd2.JPG');%读取有火焰图像
figure, subplot(121)
imshow(A)%显示A图像
subplot(122),
imshow(B)%显示B图像
%把两幅图像都分割成8x8的图像
C=cell(30,40);%把240x320图像分割成30x40个图像
for i = 1:30
    for j = 1:40
        C{i,j} = A(8*(i-1)+1:8*(i-1)+8,8*(j-1)+1:8*(j-1)+8);
    end
end
D=cell(30,40);
for i = 1:30
    for j = 1:40
        D{i,j} = B(8*(i-1)+1:8*(i-1)+8,8*(j-1)+1:8*(j-1)+8);
    end
end
%循环匹配子图像1200x1200次
for m=1:1200
    for n=1:1200
        c = normxcorr2(C{m}(:,:,1),D{n}(:,:,1));%对各块进行归一化互相关
        %figure,%创建新的图像窗口
        %surf(c),%显示两幅图像的归一化互相关
        %shading flat
        [max_c, imax] = max(abs(c(:)));%确定归一化互相关最大值及其位置
        P(m,n)=max_c;%最大值保存在P矩阵中   
    end
end


m_P_r=max(P');%对p矩阵转置求列最大值
for u=1:1200
   [row,col] = find(P(u,:)==m_P_r(u));
   zuobiao(u)=col;
end


py0 = mod(1:1200,30);%求余数（行）
for u=1:1200
    if py0(u)==0
        py0(u)=30;%把余数为0的值替换成30
    end
end
px0 = ceil((1:1200)/30); % 原始块，取整数加1（列）

py0_2=reshape(py0,30,40); % 一维转化为二维（把px0的数据转换成30x40的矩阵）
px0_2=reshape(px0,30,40); % 一维转化为二维（把px0的数据转换成30x40的矩阵）

for u=1:1200
    temp_py = mod(zuobiao(u),30);
    if temp_py==0
       temp_py=30; %把余数为0的值替换成30
    end
    temp_px = ceil(zuobiao(u)/30);
    py(u) = temp_py; % 匹配块x轴坐标
    px(u) = temp_px; % 匹配块y轴坐标
end
py_2=reshape(py,30,40); % 一维转化为二维
px_2=reshape(px,30,40);

y_p=py_2-py0_2; % y轴偏移量
x_p=px_2-px0_2; % x轴偏移量
figure
hold on
quiver(px0_2,py0_2,x_p,y_p,'-r'),hold off, axis image; %用箭头指示一个以（x_p,y_p）长度的向量


px00 = 8*px0_2-3;
py00 = 8*py0_2-3;
pxx = 8*px_2-3;
pyy = 8*py_2-3;
x_pp=px_2-px0_2; % x轴偏移量
y_pp=py_2-py0_2; % y轴偏移量