【opencv】示例-image_alignment.cpp 利用ECC 算法进行图像对齐

affine

imshow("image", target_image);
imshow("template", template_image);
imshow("warped image", warped_image);
imshow("error (black: no error)", abs(errorImage) * 255 / max_of_error);

homography

这段代码是一个利用ECC (Enhanced Correlation Coefficient) 算法进行图像对齐的示例。代码首先包含了OpenCV库的头文件，并且使用了OpenCV和标准库的命名空间。然后定义了几个函数和宏进行图像变换矩阵的操作，定义了一些用于解析命令行参数的关键字。

main 函数中，首先使用CommandLineParser解析命令行参数，之后加载输入图像和模板图像。如果只给出了输入图像，代码会生成一个随机的模板图像。如果给定了输入变换矩阵，则会使用该矩阵初始化。若变换类型为仿射或相似变换等，会初始化对应的变换矩阵。

代码接着执行findTransformECC函数来计算图像间的几何变换。计算完成后，保存变换矩阵和变换后的图像，并对结果进行显示（根据具体的verbose参数决定是否显示）。

整个过程包括了对图像的读取、几何变换矩阵的初始化与读取、使用ECC算法进行图像对齐、结果的保存与显示。通过这个示例，可以了解到图像处理中图像对齐的基本流程和方法。

/*
* 本示例展示了使用findTransformECC函数实现图像对齐的ECC算法
* 
* Demo加载一个图像（默认为fruits.jpg），并基于给定的运动类型人工创建一个模板图像。
* 当给出两幅图像时，第一幅是输入图像，第二幅定义模板图像。在后一种情况下，
* 您还可以解析warp的初始化。
* 
* 输入和输出的warp文件由原始的warp（变换）元素组成。
* 
* 作者: G. Evangelidis, 法国国家信息与自动化研究所(INRIA), Grenoble,
*       M. Asbach, 德国Fraunhofer IAIS, St. Augustin
*/// 导入opencv库的各个模块
#include <opencv2/imgcodecs.hpp> // 图像编解码
#include <opencv2/highgui.hpp>   // 图像显示
#include <opencv2/video.hpp>     // 视频处理
#include <opencv2/imgproc.hpp>   // 图像处理
#include <opencv2/core/utility.hpp> // OpenCV实用函数#include <stdio.h>  // 包含基本输入输出函数
#include <string>   // 包含字符串相关功能
#include <time.h>   // 包含处理时间相关功能
#include <iostream> // 包含输入输出流功能
#include <fstream>  // 包含文件流操作功能using namespace cv;  // 使用OpenCV命名空间
using namespace std; // 使用标准库命名空间// 向前声明函数(函数将在下面定义)
static void help(const char** argv);
static int readWarp(string iFilename, Mat& warp, int motionType);
static int saveWarp(string fileName, const Mat& warp, int motionType);
static void draw_warped_roi(Mat& image, const int width, const int height, Mat& W);// 用于向齐次变换矩阵赋值的宏定义（x和y是坐标，H是变换矩阵）
#define HOMO_VECTOR(H, x, y)\H.at<float>(0,0) = (float)(x);\H.at<float>(1,0) = (float)(y);\H.at<float>(2,0) = 1.;//  用于从齐次坐标中读取值的宏定义（X是齐次坐标，x和y是提取的坐标）
#define GET_HOMO_VALUES(X, x, y)\(x) = static_cast<float> (X.at<float>(0,0)/X.at<float>(2,0));\(y) = static_cast<float> (X.at<float>(1,0)/X.at<float>(2,0));// 定义命令行参数，用于程序运行时接收用户输入或默认值
const std::string keys ="{@inputImage    | fruits.jpg    | input image filename }" // 输入图像文件名，默认为fruits.jpg"{@templateImage |               | template image filename (optional)}" // 模板图像文件名，可选参数"{@inputWarp     |               | input warp (matrix) filename (optional)}" // 输入变换矩阵（文件名），可选参数"{n numOfIter    | 50            | ECC's iterations }" // ECC算法迭代次数，默认值为50"{e epsilon      | 0.0001        | ECC's convergence epsilon }" // ECC算法收敛阈值，默认为0.0001"{o outputWarp   | outWarp.ecc   | output warp (matrix) filename }" // 输出变换矩阵的文件名，默认为outWarp.ecc"{m motionType   | affine        | type of motion (translation, euclidean, affine, homography) }" // 变换类型，默认为仿射变换(affine)"{v verbose      | 1             | display initial and final images }" // 是否显示处理前后的图像，默认为1（显示）"{w warpedImfile | warpedECC.png | warped input image }" // 输出变形后图像的文件名，默认为warpedECC.png"{h help | | print help message }" // 显示帮助信息
;//该文件演示了 ECC 图像对齐算法的使用。 当给定一张图像时，
//模板图像是通过随机扭曲人为形成的。当给出两个图像时，可以
//通过命令行解析来初始化扭曲。如果缺少 inputWarp，则恒等变换
//会初始化算法。
// 提供程序帮助信息和命令行示例的函数
static void help(const char** argv)
{// 输出ECC图像对齐算法使用说明cout << "\nThis file demonstrates the use of the ECC image alignment algorithm. When one image"" is given, the template image is artificially formed by a random warp. When both images"" are given, the initialization of the warp by command line parsing is possible. ""If inputWarp is missing, the identity transformation initializes the algorithm. \n" << endl;// 输出仅用一个输入图像时命令行使用示例cout << "\nUsage example (one image): \n" // 使用示例（一个图像）<< argv[0]  // 程序名<< " fruits.jpg -o=outWarp.ecc "  // 输入图像，输出变换矩阵文件"-m=euclidean -e=1e-6 -N=70 -v=1 \n" << endl; // 运动类型，收敛阈值，迭代次数，是否显示结果// 输出用两个图像以及初始变换矩阵时命令行使用示例cout << "\nUsage example (two images with initialization): \n" // 使用示例（两个图像，带初始化）<< argv[0]  // 程序名<< " yourInput.png yourTemplate.png "  // 输入图像和模板图像"yourInitialWarp.ecc -o=outWarp.ecc -m=homography -e=1e-6 -N=70 -v=1 -w=yourFinalImage.png \n" << endl; // 输入变换矩阵，输出变换矩阵文件，变换类型，收敛阈值，迭代次数，是否显示结果，输出对齐后的图像文件
}// 从文件中读取变换矩阵的函数
static int readWarp(string iFilename, Mat& warp, int motionType){// 根据不同的运动类型(motionType)确定需要从文件中读取的元素数量// 如果是透视变换(Homography)，读取9个值，否则读取6个值CV_Assert(warp.type()==CV_32FC1); // 确保warp矩阵是单通道浮点型int numOfElements; // 存储元素的数量if (motionType==MOTION_HOMOGRAPHY)numOfElements=9; // 透视变换有9个参数elsenumOfElements=6; // 其它变换类型只需要6个参数int i; // 循环计数变量int ret_value; // 返回值，表示是否成功读取文件ifstream myfile(iFilename.c_str()); // 打开文件if (myfile.is_open()){ // 如果文件成功打开float* matPtr = warp.ptr<float>(0); // 获取矩阵的指针for(i=0; i<numOfElements; i++){myfile >> matPtr[i]; // 从文件读取每个值}ret_value = 1; // 读取成功，设置返回值为1}else { // 如果文件打开失败cout << "Unable to open file " << iFilename.c_str() << endl; // 输出错误信息ret_value = 0; // 设置返回值为0}return ret_value; // 返回结果
}// 将变换矩阵保存到文件中的函数
static int saveWarp(string fileName, const Mat& warp, int motionType)
{// 确保warp矩阵是单通道浮点型CV_Assert(warp.type()==CV_32FC1);const float* matPtr = warp.ptr<float>(0); // 获取矩阵的指针int ret_value; // 返回值，表示是否成功保存到文件ofstream outfile(fileName.c_str()); // 打开或创建文件来写入if( !outfile ) { // 如果文件打开失败cerr << "error in saving "<< "Couldn't open file '" << fileName.c_str() << "'!" << endl; // 输出错误信息到错误流ret_value = 0; // 设置返回值为0}else {// 如果文件成功打开, 保存warp矩阵的元素outfile << matPtr[0] << " " << matPtr[1] << " " << matPtr[2] << endl;outfile << matPtr[3] << " " << matPtr[4] << " " << matPtr[5] << endl;if (motionType==MOTION_HOMOGRAPHY){ // 如果是透视变换，需要保存额外3个参数outfile << matPtr[6] << " " << matPtr[7] << " " << matPtr[8] << endl;}ret_value = 1; // 保存成功，设置返回值为1}return ret_value; // 返回结果
}// 在图像上绘制变换后区域边界的函数
static void draw_warped_roi(Mat& image, const int width, const int height, Mat& W)
{// 定义四个角点Point2f top_left, top_right, bottom_left, bottom_right;Mat H = Mat(3, 1, CV_32F); // 定义齐次坐标向量Mat U = Mat(3, 1, CV_32F); // 存储变换后的坐标值Mat warp_mat = Mat::eye(3, 3, CV_32F); // 创建一个单位矩阵// 将变换矩阵 W 的值赋给 warp_matfor (int y = 0; y < W.rows; y++)for (int x = 0; x < W.cols; x++)warp_mat.at<float>(y,x) = W.at<float>(y,x);// 对矩形的四个角点进行变换// 左上角HOMO_VECTOR(H, 1, 1); // 将点设置为齐次坐标形式gemm(warp_mat, H, 1, 0, 0, U); // 使用矩阵乘法进行变换GET_HOMO_VALUES(U, top_left.x, top_left.y); // 获取变换后的值// 右上角HOMO_VECTOR(H, width, 1);gemm(warp_mat, H, 1, 0, 0, U);GET_HOMO_VALUES(U, top_right.x, top_right.y);// 左下角HOMO_VECTOR(H, 1, height);gemm(warp_mat, H, 1, 0, 0, U);GET_HOMO_VALUES(U, bottom_left.x, bottom_left.y);// 右下角HOMO_VECTOR(H, width, height);gemm(warp_mat, H, 1, 0, 0, U);GET_HOMO_VALUES(U, bottom_right.x, bottom_right.y);// 在图像上绘制变形后的矩形边界line(image, top_left, top_right, Scalar(255));     // 上边界line(image, top_right, bottom_right, Scalar(255)); // 右边界line(image, bottom_right, bottom_left, Scalar(255)); // 下边界line(image, bottom_left, top_left, Scalar(255));  // 左边界
}// 主函数入口
int main (const int argc, const char * argv[])
{// 解析命令行参数CommandLineParser parser(argc, argv, keys);parser.about("ECC demo"); // 关于本程序parser.printMessage(); // 打印解析的信息help(argv); // 显示帮助信息// 获取命令行参数string imgFile = parser.get<string>(0); // 输入图像文件名string tempImgFile = parser.get<string>(1); // 模板图像文件名string inWarpFile = parser.get<string>(2); // 初始变换矩阵文件名// 获取其他参数int number_of_iterations = parser.get<int>("n"); // 迭代次数double termination_eps = parser.get<double>("e"); // 精度阈值，停止条件string warpType = parser.get<string>("m"); // 变换类型int verbose = parser.get<int>("v"); // 是否显示详细信息string finalWarp = parser.get<string>("o"); // 最终变换矩阵保存文件string warpedImFile = parser.get<string>("w"); // 对齐后的图像保存文件if (!parser.check()) // 检查解析的参数是否合理{parser.printErrors();return -1; // 参数不合理则返回-1}// 确保传入的变换类型是有效的 平移、欧几里得、仿射、单应性if (!(warpType == "translation" || warpType == "euclidean"|| warpType == "affine" || warpType == "homography")){cerr << "Invalid motion transformation" << endl;return -1; // 无效变换类型，返回-1}// 根据变换类型设定变换模式int mode_temp;if (warpType == "translation")mode_temp = MOTION_TRANSLATION;else if (warpType == "euclidean")mode_temp = MOTION_EUCLIDEAN;else if (warpType == "affine")mode_temp = MOTION_AFFINE;elsemode_temp = MOTION_HOMOGRAPHY;// 读取输入图像Mat inputImage = imread(samples::findFile(imgFile), IMREAD_GRAYSCALE);if (inputImage.empty()) // 检查图像是否成功加载{cerr << "Unable to load the inputImage" <<  endl;return -1; // 加载失败，返回-1}// 初始化目标图像和模板图像Mat target_image;Mat template_image;// 如果提供了模板图像，则读取模板图像if (tempImgFile != "") {// 将输入图像复制给目标图像inputImage.copyTo(target_image);// 加载模板图像（灰度图）template_image = imread(samples::findFile(tempImgFile), IMREAD_GRAYSCALE);// 如果模板图像加载失败，则返回错误if (template_image.empty()) {cerr << "Unable to load the template image" << endl;return -1;}}else { // 如果没有指定模板图像文件名，对输入图像应用随机变换以生成模板// 对输入图像进行尺寸调整，新的尺寸为216 x 216像素resize(inputImage, target_image, Size(216, 216), 0, 0, INTER_LINEAR_EXACT);// 声明一个Mat类型变量用于存放变换矩阵Mat warpGround;// 创建一个随机数生成器，种子为当前时间戳RNG rng(getTickCount());// 声明一个double类型变量用于存放计算得到的角度double angle;// 根据mode_temp来决定如何生成变换矩阵并对图像应用变换switch (mode_temp) {case MOTION_TRANSLATION:// 生成一个随机的平移变换矩阵warpGround = (Mat_<float>(2, 3) << 1, 0, (rng.uniform(10.f, 20.f)),0, 1, (rng.uniform(10.f, 20.f)));// 应用平移变换矩阵到目标图像，获取模板图像warpAffine(target_image, template_image, warpGround,Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP);break;case MOTION_EUCLIDEAN:// 生成一个随机的欧几里得变换矩阵（含旋转和平移）angle = CV_PI / 30 + CV_PI * rng.uniform((double)-2.f, (double)2.f) / 180;warpGround = (Mat_<float>(2, 3) << cos(angle), -sin(angle), (rng.uniform(10.f, 20.f)),sin(angle), cos(angle), (rng.uniform(10.f, 20.f)));// 应用欧几里得变换矩阵到目标图像，获取模板图像warpAffine(target_image, template_image, warpGround,Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP);break;case MOTION_AFFINE:// 生成一个随机的仿射变换矩阵warpGround = (Mat_<float>(2, 3) << (1 - rng.uniform(-0.05f, 0.05f)),(rng.uniform(-0.03f, 0.03f)), (rng.uniform(10.f, 20.f)),(rng.uniform(-0.03f, 0.03f)), (1 - rng.uniform(-0.05f, 0.05f)),(rng.uniform(10.f, 20.f)));// 应用仿射变换矩阵到目标图像，获取模板图像warpAffine(target_image, template_image, warpGround,Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP);break;case MOTION_HOMOGRAPHY:// 生成一个随机的单应性变换矩阵warpGround = (Mat_<float>(3, 3) << (1 - rng.uniform(-0.05f, 0.05f)),(rng.uniform(-0.03f, 0.03f)), (rng.uniform(10.f, 20.f)),(rng.uniform(-0.03f, 0.03f)), (1 - rng.uniform(-0.05f, 0.05f)), (rng.uniform(10.f, 20.f)),(rng.uniform(0.0001f, 0.0003f)), (rng.uniform(0.0001f, 0.0003f)), 1.f);// 应用单应性变换矩阵到目标图像，获取模板图像warpPerspective(target_image, template_image, warpGround,Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP);break;}}// 根据变换类型创建适当大小的warp_matrixconst int warp_mode = mode_temp;// 初始化或加载变换矩阵Mat warp_matrix;if (warpType == "homography")// 如果是单应性变换，则warp_matrix为3x3的单位矩阵warp_matrix = Mat::eye(3, 3, CV_32F);else// 如果是其他变换，则warp_matrix为2x3的单位矩阵warp_matrix = Mat::eye(2, 3, CV_32F);// 如果提供了变换矩阵文件名，尝试读取变换矩阵if (inWarpFile != ""){int readflag = readWarp(inWarpFile, warp_matrix, warp_mode);// 如果读取失败，打印错误信息，并退出程序if ((!readflag) || warp_matrix.empty()){cerr << "-> Check warp initialization file" << endl << flush;return -1;}}else {// 如果没有提供变换矩阵文件，发出警告：假定使用单位矩阵作为初始化可能不佳，// 尤其是当图像尺寸不相同时或者形变较大时printf("\n ->Performance Warning: Identity warp ideally assumes images of ""similar size. If the deformation is strong, the identity warp may not ""be a good initialization. \n");}// 检查迭代次数是否过多if (number_of_iterations > 200)cout << "-> Warning: too many iterations " << endl;// 如果不是单应性变换，确保变换矩阵仅有两行if (warp_mode != MOTION_HOMOGRAPHY)warp_matrix.rows = 2;//开始计时const double tic_init = (double) getTickCount();// 调用findTransformECC寻找最佳变换矩阵double cc = findTransformECC(template_image, target_image, warp_matrix, warp_mode,TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, number_of_iterations, termination_eps));// 如果findTransformECC返回错误（即-1），则打印错误信息if (cc == -1){cerr << "The execution was interrupted. The correlation value is going to be minimized." << endl;cerr << "Check the warp initialization and/or the size of images." << endl << flush;}// 停止计时const double toc_final = (double)getTickCount();const double total_time = (toc_final - tic_init) / (getTickFrequency());//如果设置了verbose，打印对齐时间if (verbose) {cout << "Alignment time (" << warpType << " transformation): "<< total_time << " sec" << endl << flush;//  cout << "Final correlation: " << cc << endl << flush;}// 保存最终的变换矩阵saveWarp(finalWarp, warp_matrix, warp_mode);if (verbose) {//如果设置了verbose，打印保存变换矩阵的文件名cout << "\nThe final warp has been saved in the file: " << finalWarp << endl << flush;}// 保存最终的对齐图像Mat warped_image = Mat(template_image.rows, template_image.cols, CV_32FC1);if (warp_mode != MOTION_HOMOGRAPHY)warpAffine(target_image, warped_image, warp_matrix, warped_image.size(),INTER_LINEAR + WARP_INVERSE_MAP);elsewarpPerspective(target_image, warped_image, warp_matrix, warped_image.size(),INTER_LINEAR + WARP_INVERSE_MAP);imwrite(warpedImFile, warped_image); // 保存变形后的图像// 如果设置了verbose，显示结果图像if (verbose){cout << "The warped image has been saved in the file: " << warpedImFile << endl << flush;//创建可视化窗口namedWindow("image", WINDOW_AUTOSIZE);namedWindow("template", WINDOW_AUTOSIZE);namedWindow("warped image", WINDOW_AUTOSIZE);namedWindow("error (black: no error)", WINDOW_AUTOSIZE);//移动窗口，用于可视化布局moveWindow("image", 20, 300);moveWindow("template", 300, 300);moveWindow("warped image", 600, 300);moveWindow("error (black: no error)", 900, 300);// 绘制变换后的区域边界Mat identity_matrix = Mat::eye(3, 3, CV_32F);draw_warped_roi(target_image, template_image.cols - 2, template_image.rows - 2, warp_matrix);draw_warped_roi(template_image, template_image.cols - 2, template_image.rows - 2, identity_matrix);Mat errorImage;subtract(template_image, warped_image, errorImage);double max_of_error;minMaxLoc(errorImage, NULL, &max_of_error);// 显示图像cout << "Press any key to exit the demo (you might need to click on the images before)." << endl << flush;imshow("image", target_image);waitKey(200);imshow("template", template_image);waitKey(200);imshow("warped image", warped_image);waitKey(200);imshow("error (black: no error)", abs(errorImage) * 255 / max_of_error);waitKey(0);}// 完成程序return 0;
}

// 使用resize函数调整inputImage的大小，存入target_image中
resize(inputImage,               // 源图像target_image,             // 目标图像，输出的大小调整后的图像将被存储在这里Size(216, 216),           // 目标图像的新大小，此处指定为宽216像素，高216像素0,                        // x方向上的缩放比例，在这里缩放比例由Size参数决定，因此设置为00,                        // y方向上的缩放比例，同上设置为0INTER_LINEAR_EXACT        // 插值方式，此处使用精确线性插值算法
);

如果没有指定模板图像文件名，对输入图像应用随机变换以生成模板

warpAffine(target_image, template_image, warpGround,Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP);

// 使用增强的相关系数(ECC)算法，寻找最佳的仿射变换矩阵
double cc = findTransformECC(template_image,        // 模板图像target_image,          // 需要对齐到模板图像的对象图像warp_matrix,           // 可以指定初始估计，函数将优化这个矩阵以获得最佳变换warp_mode,             // 规定变换模型的类型，如仿射变换TermCriteria(          // 优化时的迭代终止准则，可以是最大迭代次数，变换估计的精确度，或它们的组合TermCriteria::COUNT + TermCriteria::EPS, // 使用最大迭代次数和变换估计的精确度两个条件number_of_iterations,                     // 迭代的最大次数termination_eps                            // 迭代的终止精确度)
);

// 对target_image应用透视变换，结果输出到warped_image中
warpPerspective(target_image,     // 源图像warped_image,     // 目标图像，存储变换后的图像warp_matrix,      // 透视变换矩阵，描述了变换的参数warped_image.size(), // 最终输出图像的大小INTER_LINEAR + WARP_INVERSE_MAP // 插值方式加上变换方向
);

仿射变换和透视变换详细对比

// 使用通用矩阵乘法函数gemm计算变换矩阵和透视矩阵的乘积
gemm(warp_mat,   // 第一个矩阵A，这里指的是仿射变换矩阵或相应的变换矩阵H,          // 第二个矩阵B，透视变换矩阵1,          // alpha系数，用于缩放第一个矩阵AMat(),      // 第三个矩阵C，在这个函数调用中未使用，所以传递一个空矩阵0,          // beta系数，由于矩阵C未被使用，beta系数实际不起作用U,          // 输出矩阵D，存储A和B(M1和M2)乘积的结果0           // 标志位，在此为0，默认表示正常的矩阵乘法
);