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ÖáÆ«ÒÆÁ¿