对图像作变幻，看起来是一件相对简单且有趣的事。所以，我们在选了A Look Into the Past 为主题进行实验，通过残存的光影体会照片里建筑的前世今生。
我们期望通过新旧照片特征匹配，实现图相配准(Image registration)，让旧照片融入到新照片中，达到A Look Into the Past的效果。
图像配准是使用某种算法，基于某种评估标准，将一副或多副图片（局部）最优映射到目标图片上的方法。根据不同配准方法，不同评判标准和不同图片类型，有不同类型的图像配准方法。其中，最本质的分类是：
- 基于灰度的图像配准；
- 基于特征的图像配准。
具体的图像配准算法是基于这两点的混合或者变体的算法。我们采用的是基于特征的配准，大致流程为（浏览需要用到mermaid插件）：
graph LR
A(输入图像)-->B[特征提取匹配]
B-->C[构建变换矩阵]
C-->D[新旧图像融合]
D-->E(输出)
我们采用SURF算法和KAZE算法进行特征匹配。
SURF的算法一种稳健的局部特征点检测和描述算法，其步骤如下：
graph LR
A(构建黑塞矩阵)-->B(构建尺度空间)
B-->C(特征点定位)
C-->D(特征点主方向分配)
D-->E(生成特征点描述子)
E-->F(特征点匹配)
黑塞矩阵（Hessian Matrix）是一个多元函数的二阶偏导数构成的方阵，描述了函数的局部曲率。Hessian矩阵是Surf算法的核心，构建Hessian矩阵的目的是为了生成图像稳定的边缘点（突变点），为特征提取做基础。
$$
H(f(x,y))=
\begin{bmatrix}
\frac{\partial^2f}{\partial x^2} & \frac{\partial^2f}{\partial y\partial x}\\
\frac{\partial^2f}{\partial x\partial y} & \frac{\partial^2f}{\partial y^2}
\end{bmatrix}
$$
H矩阵的判别式为
$$
det(H)=\frac{\partial^2f}{\partial x^2}\frac{\partial^2f}{\partial y^2}-(\frac{\partial^2f}{\partial y\partial x})^2
$$
当Hessian矩阵的判别式取得局部极大值时，判定当前点是比周围邻域内其他点更亮或更暗的点，由此来定位关键点的位置。
在SURF算法中，图像像素$l(x，y)$即为函数值
$$
L(x,t)=G(t)I(x,t)
$$
通过特定核间的卷积计算二阶偏导数。通过特定核间的卷积计算二阶偏导数，这样便能计算出H矩阵的三个矩阵元素$L_{xx}$,
$$
H(x,\sigma)=
\begin{bmatrix}
L_{xx}(x,\sigma) & L_{xy}(x,\sigma)\\
L_{yx}(x,\sigma) & L_{yy}(x,\sigma)
\end{bmatrix}
$$
Surf的尺度空间由O组L层组成,不同组间图像的尺寸都是一致的，但不同组间使用的盒式滤波器的模板尺寸逐渐增大，同一组间不同层间使用相同尺寸的滤波器，但是滤波器的模糊系数逐渐增大。
特征点的定位过程Surf和Sift一致，将经过Hessian矩阵处理的每个像素点与二维图像空间和尺度空间邻域内的26个点进行比较，初步定位出关键点，再经过滤除能量比较弱的关键点以及错误定位的关键点，筛选出最终的稳定的特征点。
在Surf中，采用的是统计特征点圆形邻域内的harr小波特征。 在特征点的圆形邻域内，统计60度扇形内所有点的水平、垂直harr小波特征总和，然后扇形以一定间隔进行旋转并再次统计该区域内harr小波特征值之后，最后将值最大的那个扇形的方向作为该特征点的主方向。
Surf算法中在特征点周围取一个$4\times 4$的矩形区域块，但是所取得矩形区域方向是沿着特征点的主方向。每个子区域统计25个像素的水平方向和垂直方向的haar小波特征，这里的水平和垂直方向都是相对主方向而言的。该haar小波特征为水平方向值之后、垂直方向值之后、水平方向绝对值之后以及垂直方向绝对值之和4个方向。
Surf通过计算两个特征点间的欧式距离来确定匹配度，欧氏距离越短，代表两个特征点的匹配度越好。Surf还加入了Hessian矩阵迹的判断，如果两个特征点的矩阵迹正负号相同，代表这两个特征具有相同方向上的对比度变化，如果不同，说明这两个特征点的对比度变化方向是相反的，即使欧氏距离为0，也直接予以排除。
ptsOriginal = detectSURFFeatures(original);
ptsDistorted = detectSURFFeatures(distorted);
[featuresOriginal, validPtsOriginal] = extractFeatures(original, ptsOriginal);
[featuresDistorted, validPtsDistorted] = extractFeatures(distorted, ptsDistorted);
KAZE算法是非线性扩散滤波法。具体地，非线性扩散滤波法是将图像亮度（L）在不同尺度上的变化视为某种形式的流动函数（flow function）的散度（divergence），可以通过非线性偏微分方程来描述： $$ \frac{\partial L}{\partial t}=div(c(x,y,t)\nabla L) $$ 通过设置合适的传导函数 c(x,y,t) ，可以使得扩散自适应于图像的局部结构。传导函数可以是标量、也可以是张量。时间t作为尺度参数，其值越大、则图像的表示形式越简单。Perona和Malik提出了传导函数的构造方式： $$ c(x,y,t)=g(|\nabla L_\sigma(x,y,t)|) $$ 由于非线性微分方程没有解析解，一般通过数值分析的方法进行迭代求解。传统上采用显式差分格式的求解方法只能采用小步长，收敛缓慢。KAZE中采用AOS(Additive Operator Splitting)算法对结果进行收敛。
ptsOriginal = detectKAZEFeatures(original);
ptsDistorted = detectKAZEFeatures(distorted);
这一部分的主要内容是通过匹配点对构建图像序列之间的变换矩阵。实验中通过Matlab自带的estimateGeometricTransform2D函数实现。
indexPairs = matchFeatures(featuresOriginal, featuresDistorted);
matchedOriginal = validPtsOriginal(indexPairs(:,1));
matchedDistorted = validPtsDistorted(indexPairs(:,2));
[tform, inlierIdx] = estimateGeometricTransform2D(...
matchedDistorted, matchedOriginal,'similarity','Maxdistance',1.5,'Confidence',99,'MaxNumTrials',1000);
使用imwrap函数对图像进行仿射变换，输入参数是待变换的二维图像image和一个3x3的仿射变换矩阵tform，imwrap原理：
输入图像上的点(u,v)对应输出图像上的点(x,y)，关系式为：
在计算输出图像的所有(x,y)对的时候，需要乘以A的逆矩阵然后在输入图像中插值得到输出图像。
我们的目标是实现一种从旧图像到新图像之间的渐变效果，为了达到这个效果，首先构造以旧图像中心为中心的，幅度从中心值1往边缘减小至0的函数矩阵。我们选择的是经过修正的巴特沃斯滤波器函数，定义原来图片的宽为DX，长为DY，点(x,y)的修正距离D(x,y)定义为：
巴特沃斯函数：
将$H_n (x,y)$与旧图像相乘得到输出部分1，再把$H_n (x,y)$通过相同的变换矩阵tform进行仿射变换后，取反再与新图像相乘得到输出部分2，把部分1、2相加即可得到最终结果。
f1=im2double(dist0);
t1=zeros(size(f1));
t2=t1;
LD=norm([d0-1/2,d0-1/2],2);
tx=size(t1,2);
ty=size(t1,1);
for k=1:3
for j=1:size(t1,1)
for i=1:size(t1,2)
d=norm([i/tx-0.5,j/ty-0.5],2);
t2(j,i,k)=1-1/(1+(d/LD)^p0);
t1(j,i,k)=1-t2(j,i,k);
end
end
end
dist0=im2uint8(im2double(dist0).*t1);
t1=im2uint8(t1);
t2=im2uint8(t2);
outputView = imref2d(size(original));
recovered = imwarp(dist0,tform,'OutputView',outputView);
re2=imwarp(t1,tform,'OutputView',outputView);
re2=255-re2;
orig0=im2uint8(im2double(orig0).*im2double(re2));
-
故宫
-
教堂
-
桥
.
├── code
│ ├── main.m
│ └── cv_Image_registration_.m
├── input
│ ├── new1.png
│ ├── new2.png
│ ├── new3.png
│ ├── old1.png
│ ├── old2.png
│ └── old3.png
└── output
├── result1.png
├── result2.png
└── result3.png
运行环境：
MATLAB 2021a
Computer Vision Toolbox