在前文中，我们已经有好几篇讲图像配准的文章：
3. 图像配准系列之基于B样条的FFD自由变换原理与C++实现
在以上第1、第2篇文章中，我们讲了图像配准的概念、Sift特征点的基本原理，以及基于“Sift特征点+仿射变换”的图像配准原理和C++实现。我们知道，仿射变换模型具有整体的特性，对局部形变的矫正效果并不理想。相比来说FFD变换与TPS薄板样条变换，则对局部形变具有更强的矫正能力。基于FFD变换的配准方法我们在前文已经讲过，本文我们主要讲基于“Sift特征点+TPS薄板样条变换”的配准方法，并使用C++与Opencv来将之实现。
1. Sift算法原理
在上述第1、2篇文章中已讲，此处不再重复。
2. TPS薄板样条变换原理
在上述第8篇文章中已讲，此处也不再重复。
3. 配准步骤
配准步骤的大致步骤如下图所示：
其中特征向量的匹配与剔除，我们在上方超链接的第2篇文章已讲过。
像素重采样时，由于坐标是浮点型的，会使用到插值算法，我们在之前的文章也详细说明：
4. 代码实现
好了，原理我们基本都讲过，下面直接上代码。
(1) TPS的代码
可参考开头超链接的第8篇文章，不过这里需要修正一下矩阵W的计算，调用Opencv的solve函数解矩阵方程LW=Y时，如果参数是DECOMP_LU时会出错，但是如果改为DECOMP_EIG就没问题，这里应该涉及到矩阵可逆以及求解算法的问题，有空再继续深究。
-
//W = {w0, w1, w2, ..., a1, ax, ay}
-
void cal_W(vector<Point2f> points1, vector<Point2f> points2, Mat &W)
-
{
-
int N = points1.size();
-
Mat L;
-
cal_L(points1, L);
-
-
-
Mat Y = Mat::zeros(Size(
2, N+
3), CV_32FC1);
//row:N+3, col:2
-
for(
int i =
0; i < N; i++)
-
{
-
Y.ptr<
float>(i)[
0] = points2[i].x;
-
Y.ptr<
float>(i)[
1] = points2[i].y;
-
}
-
-
-
//solve(L, Y, W, DECOMP_LU); //LW = Y,求W
-
solve(L, Y, W, DECOMP_EIG);
-
//W = L.inv()*Y;
-
}
(2) 检测特征点代码
这里有个问题需要说明一下，如果直接对整张图像进行Sift特征点检测，得到的特征点分布不均匀，然而点分布越均匀，TPS变换效果越好，因此这样直接检测整张图的特征点，最后配准效果并不理想。所以我们将图像分为多个块，分别检测各个块的特征点，然后拼接成整张图像的特征点，这样一来，特征点的分布就比较均匀了。
代码如下：
-
vector<KeyPoint> detect_sift_block(Mat img, int num, int b_row, int b_col)
-
{
-
const
int block_rsize = img.rows / b_row;
-
const
int block_csize = img.cols / b_col;
-
-
-
vector<KeyPoint> kp;
-
vector<KeyPoint> kp_tmp;
-
Ptr<SiftFeatureDetector> siftdtc = SiftFeatureDetector::create(num,
3,
0.04,
10.0,
1.6);
-
-
for (
int i =
0; i < b_row; i++)
//分块检测，确保均匀分布
-
{
-
for (
int j =
0; j < b_col; j++)
-
{
-
int xx = j*block_csize;
-
int yy = i*block_rsize;
-
Mat tmp;
-
img(Rect(xx, yy, block_csize, block_rsize)).copyTo(tmp);
-
siftdtc->detect(tmp, kp_tmp);
-
-
//将小块的点坐标转换成整张图像的点坐标
-
for (
int k =
0; k < kp_tmp.size(); k++)
-
{
-
kp_tmp[k].pt.x += xx;
-
kp_tmp[k].pt.y += yy;
-
}
-
//将多个小块的特征点拼接成整张图像的特征点
-
kp.insert(kp.end(), kp_tmp.begin(), kp_tmp.end());
-
}
-
}
-
-
-
return kp;
-
}
(3) 配准代码
使用具有局部形变的Lena图像进行配准测试，分别使用TPS变换与仿射变换作为坐标变换模型，并比较结果。
-
void Sift_Tps_test(void)
-
{
-
Mat img1 = imread(
"lena.jpg", CV_LOAD_IMAGE_GRAYSCALE);
-
Mat img2 = imread(
"lena_ffd.jpg", CV_LOAD_IMAGE_GRAYSCALE);
-
-
imshow(
"image before", img1);
-
imshow(
"image2 before", img2);
-
-
-
// SIFT - 检测关键点并在原图中绘制
-
vector<KeyPoint> kp1, kp2;
-
kp1 = detect_sift_block(img1,
50,
5,
5);
-
kp2 = detect_sift_block(img2,
50,
5,
5);
-
Mat outimg1, outimg2;
-
drawKeypoints(img1, kp1, outimg1);
-
drawKeypoints(img2, kp2, outimg2);
-
imshow(
"image1 keypoints", outimg1);
-
imshow(
"image2 keypoints", outimg2);
-
-
-
// SIFT - 特征向量提取
-
Ptr<SiftDescriptorExtractor> extractor = SiftDescriptorExtractor::create();
-
Mat descriptor1, descriptor2;
-
extractor->compute(img1, kp1, descriptor1);
-
extractor->compute(img2, kp2, descriptor2);
-
-
-
// 两张图像的特征匹配
-
Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create(
"BruteForce");
-
vector<DMatch> matches;
-
Mat img_matches;
-
matcher->match(descriptor1, descriptor2, matches);
-
-
-
//计算匹配结果中距离最大和距离最小值
-
double min_dist = matches[
0].distance, max_dist = matches[
0].distance;
-
for (
int m =
0; m < matches.size(); m++)
-
{
-
if (matches[m].distance < min_dist)
-
{
-
min_dist = matches[m].distance;
-
}
-
if (matches[m].distance > max_dist)
-
{
-
max_dist = matches[m].distance;
-
}
-
}
-
cout <<
"min dist=" << min_dist <<
endl;
-
cout <<
"max dist=" << max_dist <<
endl;
-
-
-
vector<DMatch> goodMatches;
-
for (
int i =
0; i < matches.size(); i++)
//筛选出较好的匹配点
-
{
-
if (matches[i].distance <
0.35*max_dist)
-
{
-
goodMatches.push_back(matches[i]);
-
}
-
}
-
cout <<
"The number of good matches:" << goodMatches.size() <<
endl;
-
-
-
//坐标转换为float类型
-
vector <KeyPoint> good_kp1, good_kp2;
-
for (
int i =
0; i < goodMatches.size(); i++)
-
{
-
good_kp1.push_back(kp1[goodMatches[i].queryIdx]);
-
good_kp2.push_back(kp2[goodMatches[i].trainIdx]);
-
}
-
-
-
//坐标变换
-
vector <Point2f> p01, p02;
-
for (
int i =
0; i < goodMatches.size(); i++)
-
{
-
p01.push_back(good_kp1[i].pt);
-
p02.push_back(good_kp2[i].pt);
-
}
-
-
-
vector<uchar> RANSACStatus;
//用以标记每一个匹配点的状态，等于0则为外点，等于1则为内点。
-
findFundamentalMat(p01, p02, RANSACStatus, CV_FM_RANSAC,
4.5);
//p1 p2必须为float型
-
vector<Point2f> f1_features_ok;
-
vector<Point2f> f2_features_ok;
-
for (
int i =
0; i < p01.size(); i++)
//剔除跟踪失败点
-
{
-
if (RANSACStatus[i])
-
{
-
f1_features_ok.push_back(p01[i]);
//基准图特征点
-
f2_features_ok.push_back(p02[i]);
//流动图特征点
-
}
-
}
-
-
-
#if 1 //这里使用TPS变换
-
Mat Tx, Ty;
-
Tps_TxTy(f1_features_ok, f2_features_ok, img2.rows, img2.cols, Tx, Ty);
-
Mat tps_out;
-
remap(img2, tps_out, Tx, Ty, INTER_CUBIC);
-
-
-
imshow(
"img2-img1",
abs(img2-img1));
-
imshow(
"tps_out-img1",
abs(tps_out - img1));
-
imshow(
"tps_out", tps_out);
-
#else //这里使用仿射变换
-
Mat T = estimateRigidTransform(f2_features_ok, f1_features_ok,
true);
//false表示刚性变换
-
-
-
Mat affine_out;
-
warpAffine(img2, affine_out, T, img2.size(), INTER_CUBIC);
-
-
-
imshow(
"img2-img1",
abs(img2 - img1));
-
imshow(
"affine_out-img1",
abs(affine_out - img1));
-
imshow(
"affine_out", affine_out);
-
#endif
-
-
-
cv::waitKey(
0);
-
}
5. 配准结果
(1) 首先是使用“Sift+仿射变换”的结果：
参考图像
浮动图像
参考图像的Sift特征点
浮动图像的Sift特征点
配准图像
参考图像与浮动图像的差值图
参考图像与配准图像的差值图
(2) 其次是使用“Sift+TPS变换”的结果：
配准图像
参考图像与浮动图像的差值图
参考图像与配准图像的差值图
由以上结果可知，仿射变换模型并不能很好的矫正局部形变，相比来说，TPS变换却可以较好地矫正局部形变。不过TPS计算更复杂更耗时，因此在实际操作中，我们需要根据图像地特点来选择合适的坐标变换模型，使配准效果与配准耗时都满足应用场景要求。
欢迎扫码关注以下微信公众号，接下来会不定时更新更加精彩的内容噢～
转载：https://blog.csdn.net/shandianfengfan/article/details/114386128