MATLAB学习笔记——自相关/互相关函数xcorr使用介绍

原文链接：https://blog.csdn.net/lvsehaiyang1993/article/details/82913443

转载请附上原文出处链接

侵删

Matlab提供了计算互相关和自相关的函数xcorr函数
1.使用方法

c = xcorr(x,y)
c = xcorr(x)
c = xcorr(x,y,‘option’)
c = xcorr(x,‘option’)
c = xcorr(x,y,maxlags)
c = xcorr(x,maxlags)
c =
xcorr(x,y,maxlags,‘option’)
c =
xcorr(x,maxlags,‘option’)
[c,lags] = xcorr(…)
其中option为：
"biased"为有偏的互相关函数估计；
"unbiased"为无偏的互相关函数估计；
"coeff"为0延时的正规化序列的自相关计算；
"none"为原始的互相关计算
2.特别的：
c=
xcorr(x,‘option’)特指以上某个选项的自相关估计。
c = xcorr(x,y,maxlags)
返回一个延迟范围在[-maxlags,maxlags]的互相关函数序列，输出c的程度为2maxlags+1.
c = xcorr(x,maxlags)
返回一个延迟范围在[-maxlags,maxlags]的自相关函数序列，输出c的程度为2maxlags+1.
c = xcorr(x,y,maxlags,‘option’)
同时指定maxlags和option的互相关计算.
c = xcorr(x,maxlags,‘option’)
同时指定maxlags和option的自相关计算.
3.例子
3.1计算自相关

dt=.1;
t=[0:dt:100];
x=cos(t);
[a,b]=xcorr(x,‘coeff’);
plot(b*dt,a)

3.2计算互相关

dt=.1;
t=[0:dt:1000];
x=cos(t);
y=sin(t);
[a,b]=xcorr(x,y,200,‘coeff’);
plot(b*dt,a)

————————————————
原文链接：https://blog.csdn.net/lvsehaiyang1993/article/details/82913443

转载请附上原文出处链接

零均值归一化互相关函数相对于别的一些相关函数，精度虽然没有多大提升，但是它具有较强的抗干扰性和鲁棒性。所以，在数字图像相关法（DIC）的图像匹配过程中，一般都是使用零均值归一化互相关函数，下面是其公式：

在进行数字图像相关匹配时，通常分为两部分，首先是进行整像素的定位，然后是亚像素的精确计算。在进行亚像素的精确求解时，需要以整像素的定位结果为基础，因此，一个准确的整像素匹配坐标对于后续的计算来说至关重要。数字图像相关方法发展至今，已有许多学者提出了多种整像素匹配方法，例如初期的逐点法，以及后期的爬山法、粒子群法等智能算法。虽然逐点法相对于之后的一些智能算法在速度上要落后很多，但是由于它是一种全局搜索方法，因此在准确性上要大大优于其它的方法。本节主要介绍了一种快速 ZNCC 算法，不仅保证了传统逐点算法的准确性，还大大提高了计算速度（摘抄自《基于双目视觉的DIC测量技术研究--肖颖》第三章3.2整像素相关匹配算法）。快速ZNCC算法原文请查看我的上传资源。

下面是MATLAB代码实现：

这是主脚本：

%% 快速零均值归一化互相关函数
%结果为1相关性最好，值越小相关性最差
%计算两个大小相同的灰度矩阵的相关系数
%%
tic
clear;
clc;
[filename1,pathname1,~]=uigetfile({'*.jpg';'*.jpeg';'*.bmp'},'选择一张图片');
orignf=imread([pathname1 filename1]);    %读取参考子区色度信息
f=rgb2gray(orignf);                    %将照片转化为灰度照片
f=double(f);
[r1,c1]=size(f);
if rem(r1,2)==0
    f=imresize(f,[r1-1,c1]);
elseif rem(c1,2)==0
    f=imresize(f,[r1,c1-1]);
end
[r1,c1]=size(f);
[filename2,pathname2,~]=uigetfile({'*.jpg';'*.jpeg';'*.bmp'},'选择一张图片');
orignG=imread([pathname2 filename2]);    %读取目标图像照片色度信息
G=rgb2gray(orignG);                    %读取照片色度信息
G=double(G);
[r2,c2]=size(G);
%%
F=sum(sum(f))/numel(f);           %参考子区照片灰度值平均值
den1=0;                           %分母左部分
for i=1:r1
    for j=1:c1
        den1=den1+(f(i,j)-F)^2;
    end
end
den1=sqrt(den1);
Subf=0;
for i=1:r1
    for j=1:c1
        Subf=Subf+f(i,j);
    end
end
m=r1;
n=c1;
Xmin=(m-1)/2+2;       %子区中心可选范围
Xmax=r2-(m-1)/2;
Ymin=(n-1)/2+2;
Ymax=c2-(n-1)/2;
Sg=SumTableg(G,r2,c2);   %累加和表
Sgg=SumTablegg(G,r2,c2);
COE=zeros(Xmax-Xmin+1,Ymax-Ymin+1);
num=0;
den2=0;
for x=Xmin:Xmax
    for y=Ymin:Ymax
        g=G(x-(m-1)/2:x+(m-1)/2,y-(n-1)/2:y+(n-1)/2);    %目标子区
        Sfg=SumTablefg(f,g,r1,c1);
        Meang=sum(sum(g))/numel(g);           %目标子区平均灰度值
        Subg=Sg(x+(m-1)/2,y+(n-1)/2)-Sg(x-(m-1)/2-1,y+(n-1)/2)...
            -Sg(x+(m-1)/2,y-(n-1)/2-1)+Sg(x-(m-1)/2-1,y-(n-1)/2-1);
        Subgg=Sgg(x+(m-1)/2,y+(n-1)/2)-Sgg(x-(m-1)/2-1,y+(n-1)/2)...
            -Sgg(x+(m-1)/2,y-(n-1)/2-1)+Sgg(x-(m-1)/2-1,y-(n-1)/2-1);
        Subfg=Sfg(end,end);
        num=Subfg-Meang*Subf-F*Subg+r1*c1*F*Meang;
        den2=Subgg-2*Meang*Subg+Meang*Meang*r1*c1;
        den2=sqrt(den2);
        C=num/(den1*den2);
        COE(x-Xmin+1,y-Ymin+1)=C;
    end
end
%% 
[M,~]=max(COE);
[~,J]=max(M);
y=Ymin+J;
[M,I]=max(COE(:,J));
x=Xmin+I;
fprintf('最大相关系数是：%.6f\n',M)
L=2;   %画框宽度，单位像素
orignGCopy(:,:,1)=orignG(x-(m-1)/2+L:x+(m-1)/2-L,y-(n-1)/2+L:y+(n-1)/2-L,1);
orignGCopy(:,:,2)=orignG(x-(m-1)/2+L:x+(m-1)/2-L,y-(n-1)/2+L:y+(n-1)/2-L,2);
orignGCopy(:,:,3)=orignG(x-(m-1)/2+L:x+(m-1)/2-L,y-(n-1)/2+L:y+(n-1)/2-L,3);
%% 
orignG(x-(m-1)/2:x+(m-1)/2,y-(n-1)/2:y+(n-1)/2,1)=255;
orignG(x-(m-1)/2:x+(m-1)/2,y-(n-1)/2:y+(n-1)/2,2)=0;
orignG(x-(m-1)/2:x+(m-1)/2,y-(n-1)/2:y+(n-1)/2,3)=0;
orignG(x-(m-1)/2+L:x+(m-1)/2-L,y-(n-1)/2+L:y+(n-1)/2-L,1)=orignGCopy(:,:,1);
orignG(x-(m-1)/2+L:x+(m-1)/2-L,y-(n-1)/2+L:y+(n-1)/2-L,2)=orignGCopy(:,:,2);
orignG(x-(m-1)/2+L:x+(m-1)/2-L,y-(n-1)/2+L:y+(n-1)/2-L,3)=orignGCopy(:,:,3);

image(orignG)
toc

下面的这些都是主脚本所需函数：

function GraySum=SumTableg(P,m,n)
S=zeros(m,n);
for x=1:m
    for y=1:n
        if x==1
            for j=1:y
                S(x,y)=S(x,y)+P(x,j);
            end
        elseif y==1&&x~=1
            for i=1:x
                S(x,y)=S(x,y)+P(i,y);
            end
        else
            S(x,y)=P(x,y)+S(x-1,y)+S(x,y-1)-S(x-1,y-1);
        end
    end
end
GraySum=S;
end

function GraySum=SumTablegg(P,m,n)
S=zeros(m,n);
for x=1:m
    for y=1:n
        if x==1
            for j=1:y
                S(x,y)=S(x,y)+P(x,j).^2;
            end
        elseif y==1&&x~=1
            for i=1:x
                S(x,y)=S(x,y)+P(i,y).^2;
            end
        else
            S(x,y)=P(x,y).^2+S(x-1,y)+S(x,y-1)-S(x-1,y-1);
        end
    end
end
GraySum=S;
end

function GraySum=SumTablefg(F,G,m,n)
S=zeros(m,n);
for x=1:m
    for y=1:n
        if x==1
            for j=1:y
                S(x,y)=S(x,y)+F(x,j)*G(x,j);
            end
        elseif y==1&&x~=1
            for i=1:x
                S(x,y)=S(x,y)+F(i,y)*G(i,y);
            end
        else
            S(x,y)=F(x,y)*G(x,y)+S(x-1,y)+S(x,y-1)-S(x-1,y-1);
        end
    end
end
GraySum=S;
end

在运行程序时，选择的第一张图片是参考子区（像素少的），第二张是目标图像（像素多的），千万不要错了！！！

下面是两张图片：

                                       图一                                                 图二