一、基于灰度的模板匹配

参考文章：https://blog.csdn.net/hujingshuang/article/details/47759579
【加*表示可适用于多模图像配准】

1.1 MAD算法

平均绝对差算法（Mean Absolute Differences，简称MAD算法）。它是Leese在1971年提出的一种匹配算法。算法的思想简单，具有较高的匹配精度。

算法思路： 选模板$R_{m\times n}$，从搜索图$S_{M\times N}$左上角像素开始在$(M-m,N-n)$的范围内逐行逐列进行搜索，寻找相似度最高的子图$s_{m\times n}$。

相似度定义为：
$$D(i,j)=\frac{1}{M\times N}\sum _{s=1}^M\sum _{t=1}^N|S(i+s-1,j+t-1)-T(s,t)|，1\le i\le m-M+1,1\le j\le n-N+1$$

1.2 SAD算法(与MAD类似）

绝对误差和算法（Sum of Absolute Differences，简称SAD算法）。实际上，SAD算法与MAD算法思想几乎是完全一致，只是其相似度测量公式有一点改动（计算的是子图与模板图的L1距离）

相似度定义为：
$$D(i,j)=\sum _{s=1}^M\sum _{t=1}^N|S(i+s-1,j+t-1)-T(s,t){|}_1$$

matlab代码实现：

src=imread('lena.jpg');
[M,N,d]=size(src);
if d==3src=rgb2gray(src);
end
mask=imread('lena2.jpg');
[m,n,d]=size(mask);
if d==3mask=rgb2gray(mask);
enddst=zeros(M-m,N-n); %//保存相似度
for i=1:M-m         %//子图选取，每次滑动一个像素for j=1:N-ntemp=src(i:i+m-1,j:j+n-1);%//当前子图dst(i,j)=dst(i,j)+sum(sum(abs(temp-mask)));end
end
abs_min=min(min(dst));
[x,y]=find(dst==abs_min);%//寻找相似度最高（即D(i,j)最小）的位置
figure;
imshow(mask);title('模板');
figure;
imshow(src);
hold on;
rectangle('position',[y,x,n,m],'edgecolor','r');
hold off;title('搜索图');

1.3 SSD算法

平均误差平方和算法（Mean Square Differences，简称MSD算法），也称均方差算法。

相似度定义为：
$$D(i,j)=\sum _{s=1}^M\sum _{t=1}^N[S(i+s-1,j+t-1)-T(s,t){]}^2$$其中$T(s,t)$表示子图的均值。

1.4 MSD算法（与SSD类似）

平均误差平方和算法（Mean Square Differences，简称MSD算法），也称均方差算法。

相似度定义为：
$$D(i,j)=\frac{1}{M\times N}\sum _{s=1}^M\sum _{t=1}^N[S(i+s-1,j+t-1)-T(s,t){]}^2$$

1.5 NCC算法

归一化积相关算法（Normalized Cross Correlation，简称NCC算法），与上面算法相似，依然是利用子图与模板图的灰度，通过归一化的相关性度量公式来计算二者之间的匹配程度。

相似度定义为：
$$R(i,j)=\frac{{\sum }_{s=1}^M{\sum }_{t=1}^N\left|S^{i,j}(s,t)-E\left(S^{i,j}\right)\right|\cdot |T(s,t)-E(T))|}{\sqrt{{\sum }_{s=1}^M{\sum }_{t=1}^N{\left[S^{i,j}(s,t)-E\left(S^{i,j}\right)\right]}^2\cdot {\sum }_{s=1}^M{\sum }_{t=1}^N[T(s,t)-E(T){]}^2}}$$其中$S^{i,j}(s,t)$和$E(T)$分别表示子图和模板的均值。

matlab代码：

function [ncc] = NCC(m1,m2) 
% 用于计算两幅图像的NCC，m1是基准[M,N]=size(m1);[m,n]=size(m2);ES=sum(sum(m1))/(M*N);ET=sum(sum(m2))/(m*n);A = abs(double(m1)-ES);B = abs(double(m2)-ET);AB = sum(sum(A.*B));C = A.^2;D = B.^2;CD = sum(sum(C)).*sum(sum(D));ncc = AB/sqrt(CD);
end

1.6 SSDA算法

序贯相似性检测算法（Sequential Similiarity Detection Algorithm，简称SSDA算法），它是由Barnea和Sliverman于1972年，在文章《A class of algorithms for fast digital image registration》中提出的一种匹配算法，是对传统模板匹配算法的改进，比MAD算法快几十到几百倍。

算法思路： 分别从子图和模板中随机选择一些对应像素点进行相似度比较，从而加快检测的速度。

算法步骤：

1. 设置阈值$T_h$（阈值越大越准确，但是越慢），遍历子图，求模板图和子图的图像均值
   $\overline{S_{l,j}}=E\left(S_{i,j}\right)=\frac{1}{M\times N}{\sum }_{s=1}^M{\sum }_{t=1}^N{S}_{i,j}(s,t)$ ， $\overline{T}=E(T)=\frac{1}{M\times N}{\sum }_{s=1}^M{\sum }_{t=1}^{N}T(s,t)$

2. 随机取不重复的点，求模板图和子图中该点的值与其自身图像均值的绝对误差：$\epsilon (i,j,s,t)=\left|S_{i,j}(s,t)-\overline{S_{l,j}}-(T(s,t)-\overline{T})\right|$

3. 累加绝对误差，直到超过阈值$T_h$，统计累加次数

4. 累加次数最多的子图为最佳匹配的图像【若最大的累加次数有多个（一般不存在），则取累加误差最小的作为匹配图像】

python代码实现：
【参考文章】https://www.jianshu.com/p/ecf44fa8a8ba

import time
import random
import cv2
import numpy as npdef Image_Mean_Value(Image):res = np.sum((Image.astype("float"))) / float(Image.shape[0] * Image.shape[1])return resdef Image_SSDA(search_img, example_img):# 确定子图的范围M = search_img.shape[0]m = example_img.shape[0]N = search_img.shape[1]n = example_img.shape[1]Range_x = M - m - 1Range_y = N - n - 1search_img = cv2.cvtColor(search_img, cv2.COLOR_RGB2GRAY)example_img = cv2.cvtColor(example_img, cv2.COLOR_RGB2GRAY)# 求模板图的均值example_MV = Image_Mean_Value(example_img)R_count_MAX = 0Best_x = 0Best_y = 0# 遍历所有可能的子图for i in range(Range_x):for j in range(Range_y):R_count = 0Absolute_error_value = 0# 从搜索图中截取子图subgraph_img = search_img[i:i+m, j:j+n]# 求子图的均值search_MV = Image_Mean_Value(subgraph_img)# 判断是否超过阈值while Absolute_error_value < SSDA_Th:R_count += 1# 取随机的像素点x = random.randint(0, m - 1)y = random.randint(0, n - 1)pv_s = subgraph_img[x][y]pv_e = example_img[x][y]# 计算绝对误差值并累加Absolute_error_value += abs((pv_e - example_MV) - (pv_s - search_MV)) # 如果子图和模板完全一样，绝对误差累加一直为0会陷入死循环if R_count>m*n:break  # 代码改进：规定累加次数不超过模板的像素数# 将次数最多的子图为匹配图像if R_count >= R_count_MAX:R_count_MAX = R_countBest_x = iBest_y = j# 返回最匹配图像的坐标return Best_x, Best_yif __name__ == '__main__':# 原图路径srcImg_path = "lena.jpg"# 搜索图像路径searchImg_path = "lena2.jpg"# SSDA算法阈值SSDA_Th = 10src_img = cv2.imread(srcImg_path)search_img = cv2.imread(searchImg_path)start = time.perf_counter()Best_x, Best_y = Image_SSDA(src_img, search_img)end = time.perf_counter()print("time:", end - start)cv2.rectangle(src_img, (Best_y, Best_x), (Best_y + search_img.shape[1], Best_x + search_img.shape[0]), (0, 0, 255), 1)cv2.imshow("src_img", src_img)cv2.imshow("search_img", search_img)cv2.waitKey(0)

1.7 SATD算法

hadamard变换算法（Sum of Absolute Transformed Difference，简称SATD算法），它是经hadamard变换再对绝对值求和算法。

算法思路： 将模板与子图做差后得到的矩阵Q，再对矩阵Q求其hadamard变换（左右同时乘以H，即HQH），对变换都得矩阵求其元素的绝对值之和即SATD值，作为相似度的判别依据，SATD值最小的子图，便是最佳匹配。

matlab代码实现：

src=double(rgb2gray(imread('lena.jpg')));%//长宽相等的
mask=double(rgb2gray(imread('lena3.jpg')));%//长宽相等的
M=size(src,1);%//搜索图大小
N=size(mask,1);%//模板大小
hdm_matrix=hadamard(N);%//hadamard变换矩阵
hdm=zeros(M-N,M-N);%//保存SATD值
for i=1:M-Nfor j=1:M-Ntemp=(src(i:i+N-1,j:j+N-1)-mask)/256;sw=(hdm_matrix*temp*hdm_matrix)/256;hdm(i,j)=sum(sum(abs(sw)));end
end
min_hdm=min(min(hdm));
[x,y]=find(hdm==min_hdm);
figure;imshow(uint8(mask));
title('模板');
figure;imshow(uint8(src));hold on;
rectangle('position',[y,x,N,N],'edgecolor','r');
title('搜索结果');hold off;

存在的问题： 模板由于要进行hadamard变换，所以需要选取$2^n\ast 2^n$大小的像素块进行匹配。

1.8 局部灰度值编码算法

【参考文章】https://hujingshuang.blog.csdn.net/article/details/47803791
算法思路： 通过对灰度值编码来进行粗匹配，再用相位相关法进行精匹配。

粗匹配： 略

精匹配： 使用相位相关法，根据二维傅里叶变换的性质：空域上的平移等价于频域相位的平移。两幅图的平移矢量可以通过他们的互功率谱的相位直接计算。具体计算过程：

1. 假设额图像$f_1(x,y)$和$f_2(x,y)$之间的位移矢量为$(x_0,y_0)$，即：${\mathrm{f}}_1(\mathrm{x},\mathrm{y})={\mathrm{f}}_2\left(\mathrm{x}-{\mathrm{x}}_0,\mathrm{y}-{\mathrm{y}}_0\right)$
2. 分别对$f_1(x,y)$和$f_2(x,y)$进行傅里叶变换，在频域中$F_1(u,v)$和$F_2(u,v)$的关系为：$F_1(u,v)=F_2(u,v)e^{-j2\pi \left(u{x}_0+v{y}_0\right)}$
3. 求它们的归一化互功率谱为：$\frac{F_1(u,v)F_2^{\ast }(u,v)}{|F_1(u,v)F_2^{\ast }(u,v)|}=e^{-j2\pi \left(u{x}_0+v{y}_0\right)}$，其中 $F^{\ast }$为$F$的复数共轭
4. 对上述归一化互功率谱等式左右两边进行二维傅里叶逆变换，右边逆变换可得到在空间位置$(x_0,y_0)$处形成的一个脉冲函数。所以对左边式子进行逆变换，再寻找脉冲位置（即峰值），峰值的位置$(x_0,y_0)$就是偏移矢量。

matlab代码实现精匹配：

img=rgb2gray(imread('ena.jpg'));%//读取搜索图像
dl=64;%//模板大小64x64
x0=200;y0=200;%//定义起点
dx0=20;dy0=30;%//设置偏移量（不超过模板尺寸）
f1=img(x0:x0+dl-1,y0:y0+dl-1);%//截取一个图当做模板
f2=img(x0+dx0:x0+dx0+dl-1,y0+dy0:y0+dy0+dl-1);%//截取一个图当做粗匹配找到的子图
[m,n]=size(f1);%//m行n列
F1=fft2(f1,m,n);%//傅里叶变换
F2=fft2(f2,m,n);
%//FC1=conj(F1);%//共轭
FC2=conj(F2);%//共轭
delta=abs(ifft2((F1.*FC2)./abs(F1.*FC2))); %//归一化互功率谱
figure;surf(1:dl,1:dl,delta);shading interp; %//绘制峰图
max_delat=max(max(delta)); %//峰值
[dx,dy]=find(delta==max_delat); %//dx,dy是矩阵坐标，从1开始
dx=dx-1;dy=dy-1; %//转化到从0开始

1.9 *PIU算法

划分强度一致法（Partitioned Intensity Uniformity，PIU），是针对多模图像匹配的一种算法，例如医学图像中，同一物体不同的模态下灰度值是不同的，此时传统相似度衡量的算法就失效了。PIU定义为：
$$\mathrm{PIU}(R,F)=\sum _r\frac{n_r{\sigma }_{F_{\tau }}(r)}{N{\mu }_{F_{\tau }}(r)}+\sum _f\frac{n_f{\sigma }_R(f)}{N{\mu }_R(f)}$$其中$n_r,n_f$分别表示模板和子图中灰度为$r,f$的像素数，$F_r$表示某个子图，$N$表示模板中的总像素数，另外：
$$\begin{gathered} {\mu }_{F_{\tau }}(r)=\frac{1}{n_r}\sum _{{\Omega }_r}{F}_{\tau }\left(x_R\right) \\ {\mu }_R(f)=\frac{1}{n_f}\sum _{{\Omega }_f}R\left(x_{F_{\tau }}\right) \\ {\sigma }_{F_{\tau }}(r)=\frac{1}{n_r}\sum _{{\Omega }_r}{\left[F_{\tau }\left(x_R\right)-{\mu }_{F_{\tau }}(r)\right]}^2 \\ {\sigma }_R(r)=\frac{1}{n_f}\sum _{{\Omega }_f}{\left[R\left(x_{F_{\tau }}\right)-{\mu }_R(f)\right]}^2 \end{gathered}$$其中${\sum }_{{\Omega }_r}{F}_{\tau }\left(x_R\right)$表示模板$R$中灰度值为$r$的像素在子图中对应位置上像素灰度值之和。其他式子也是同样的道理。

缺点： 复杂度很高

matlab代码实现：

img=rgb2gray(imread('C:\Users\nice\Pictures\lena.jpg'));
[img_rows,img_clos]=size(img);%//搜索图尺寸
x0=256;y0=256;len=64;
R=img(x0:x0+len-1,y0:y0+len-1);%//取一部分作为模板
R=histeq(R); %对模板图像进行均衡化，改变灰度分布
figure;imshow(R);title('模板');
[rows,clos]=size(R);%//模板尺寸
N=rows*clos;
uFr=zeros(256,1);uRf=zeros(256,1);
dFr=zeros(256,1);dRf=zeros(256,1);
a=zeros(256,1);b=zeros(256,1);
piu=zeros(img_rows-len,img_clos-len);
for i=1:img_rows-lenfor j=1:1:img_clos-len %//原文代码中将步长设置为5跑不通，改为正常步长就好S=img(i:i+len-1,j:j+len-1);%子图for r=0:255pos=find(R==r);nr=size(pos,1)+eps;%//参考图像中灰度值为r的像素个数,+eps是为了防止nf为零后面做除数出错value1=S(pos);uFr(r+1,1)=sum(value1)/nr;% //R中像素为r的位置处，对应S上的像素点均值t1=double(S(pos))-uFr(r+1,1);dFr(r+1,1)=sum(t1.^2)/nr;pos=find(S==r);nf=size(pos,1)+eps;%//参考图像中灰度值为r的像素个数value2=R(pos);uRf(r+1,1)=sum(value2)/nf;t2=double(R(pos))-uRf(r+1,1);dRf(r+1,1)=sum(t2.^2)/nf;a(r+1,1)=(nr*dFr(r+1,1))/(N*uFr(r+1,1)+eps);b(r+1,1)=(nf*dRf(r+1,1))/(N*uRf(r+1,1)+eps);endpiu(i,j)=sum(a)+sum(b);end
end
piu_min=min(min(piu));
[y,x]=find(piu==piu_min);
x=x-1;y=y-1;
figure;imshow(img);hold on;
rectangle('position',[y,x,len-1,len-1],'edgecolor','r');
title('搜索结果');hold off;

2.0 *基于互信息的算法（MI、NMI、ECC）

基于互信息的医学图像配准方法（Mutual Information，MI）被认为是最好的配准方法之一。

两幅图的互信息是通过它们的熵以及联合熵，来反映它们之间信息的相互包含程度。对于图像X、Y来说，其互信息表示为：
$$MI(X,Y)=\sum _{x\in X}\sum _{y\in Y}p(x,y)\log \frac{p(x,y)}{p(x)p(y)}$$熵定义为：
$$H(X)=-\sum _{x\in X}p\left(x\right)\log p\left(x\right)$$联合熵定义为：
$$H(X,Y)=-\sum _{x,y}p(x,y)\log p(x,y)$$其中$p(x,y)$表示联合分布，$p(x),p(y)$是边缘分布，则互信息公式可以简单表示为：
$$\begin{array}{rl}MI(X,Y)& =H(X)-H(X\mid Y)\\ & =H(X)+H(Y)-H(X,Y)\end{array}$$寻找模板与各子图之间互信息（MI）的最大值，即为配准图像。

互信息改进算法： 互信息对两幅图像之间的重叠区域比较敏感，如果两幅图像的重叠区太小，互信息就会很小，配准精度随之降低。改进后提出了归一化互信息（Normalization Mutual Information，NMI）、熵相关系数（Entropy Corrleation Coefficient，ECC），表达式如下：
$$\begin{array}{rl}NM(X,Y)& =\frac{H(X)+H(Y)}{H(X,Y)}\\ ECC(X,Y)& =\frac{2MI(X,Y)}{H(X)+H(Y)}\end{array}$$

matlab代码实现：

img=rgb2gray(imread('lena.jpg'));%//搜索图
[M,N]=size(img);
%//--------------------------------------------------------------------------
x0=256;y0=256;%选取模板起点的位置
dx=64;dy=64;%//模板、子图尺寸
img1=img(x0:x0+dx-1,y0:y0+dy-1);%//模板
img1=histeq(img1);
ET=entropy(img1);%//模板熵
%//--------------------------------------------------------------------------
%//联合熵
[m,n]=size(img1);%//模板尺寸
MI=zeros(M-dx,N-dy);%//互信息
NMI=zeros(M-dx,N-dy);%//归一化互信息
ECC=zeros(M-dx,N-dy);%//熵相关系数
for i=1:M-dxfor j=1:N-dyimg2=img(i:i+dx-1,j:j+dy-1);%//子图ES=entropy(img2);%//模板熵histq=zeros(256,256);%//联合直方图，清空%//联合直方图for s=1:mfor t=1:nx=img1(s,t)+1;y=img2(s,t)+1;%//灰度<—>坐标histq(x,y)=histq(x,y)+1;endendp=histq./sum(sum(histq));%//联合概率密度EST=-sum(sum(p.*log(p+eps)));%//联合熵（越小说明相似度越高）MI(i,j)=ES+ET-EST;%//MI互信息越大，说明相互包含的信息多，即越匹配NMI(i,j)=(ES+ET)/EST;%//NMI,越大越匹配ECC(i,j)=2*MI(i,j)/(ES+ET);%//ECC，越大越匹配end
end
%//--------------------------------------------------------------------------
mi_max=max(max(MI));
nmi_max=max(max(NMI));
ncc_max=max(max(ECC));
[xt1,yt1]=find(MI==mi_max);
[xt2,yt2]=find(NMI==nmi_max);
[xt3,yt3]=find(ECC==ncc_max);
src=img1;
dst1=img(xt1:xt1+dx-1,yt1:yt1+dx-1);
dst2=img(xt2:xt2+dx-1,yt2:yt2+dx-1);
dst3=img(xt3:xt3+dx-1,yt3:yt3+dx-1);
figure;imshow(src);title('模板');
figure;imshow(img);hold on;rectangle('position',[yt1,xt1,n-1,m-1],'edgecolor','r');title('MI配准图');hold off;
figure;imshow(img);hold on;rectangle('position',[yt2,xt2,n-1,m-1],'edgecolor','r');title('NMI配准图');hold off;
figure;imshow(img);hold on;rectangle('position',[yt3,xt3,n-1,m-1],'edgecolor','r');title('NCC配准图');hold off;