function jieguo= PC_O(cut,tar)
%UNTITLED2 此处显示有关此函数的摘要
% 此处显示详细说明
%PC_O(phase correlation operated in the frequency domain using FFT on orientation images)
esnr=eSnr(cut);%估计图像信噪比
Xsobel=[-1,0,1;-2,0,2;-1,0,1;];%利用sobel算子求梯度（差分）
Ysobel=Xsobel';
XOcut=conv2(cut,Xsobel,'same');%图像卷积
YOcut=conv2(cut,Ysobel,'same');%图像卷积
Ocut=sign(XOcut+YOcut*i);%计算orientation image
XOtar=conv2(tar,Xsobel,'same');%图像卷积
YOtar=conv2(tar,Ysobel,'same');%图像卷积
Otar=sign(XOtar+YOtar*i);%计算orientation image
[trows,tcols]=size(tar);
fftcut=fft2(Ocut,trows,tcols);%二维离散傅里叶变换,变成mrows*ncols模板,原来图像的部分填0，原来图像在左上角
ffttar=fft2(Otar);%二维离散傅里叶变换
cfftcut=conj(fftcut);%共轭计算
ff=cfftcut.*ffttar;%功率谱
fff=abs(ff);%卧槽 fff团
c=ff./fff;%相位谱
cc=ifft2(ff);%反傅里叶变换
[max_c, imax] = max(cc(:));%找到最大值，（：）表示将矩阵以一列输出，max_c是最大值 imax是下标
[xpeak, ypeak] = ind2sub(size(cc),imax(1));%返回索引值imax(1)所在的行列号
if xpeak-2<=0 || xpeak+2>size(cc,1) || ypeak-2<=0 || ypeak+2>size(cc,2)
jieguo=[0,0,0,0;];%如果峰值过于靠近边缘，则认为是误匹配
else
% fit=surface_fitting(cc(xpeak-2:xpeak+2,ypeak-2:ypeak+2));%亚像素匹配
jieguo= [xpeak,ypeak,max_c,esnr];
%上面加1是因为两个矩阵卷积的结果是两个矩阵的行列大小之和减1,例如6*6的矩阵和6*6的矩阵卷积就是11*11的矩阵 那么
%为了得到峰值点在原矩阵上的位置，首先要减去搜索窗的大小,但减去之后还要加上1,因为上述的减1的原因。
end
end
