OpenCV + CPP 系列（廿三）像素重映射 与 图像扭曲(MLS)

一、重映射简介

重映射，就是把一幅图像中某位置的像素放置到另一图像指定位置的过程。即：
$$\mathrm{d}\mathrm{s}\mathrm{t}(x,y)=\mathrm{s}\mathrm{r}\mathrm{c}(x\_map(x,y),y\_map(x,y))$$

在重映射过程中，图像的大小也可以同时发生改变。此时像素与像素之间的关系就不是一一对应关系，因此在重映射过程中，可能会涉及到像素值的插值计算。点击查看边缘处理项

Remap(

InputArray src,　　　　　　　输入图像（灰度图或真彩图均可）
OutputArray dst,　　　 　　　输出图像（要求大小和xmap,ymap相同，通道数目及数据类型和src相同）
InputArray map1,　　　　　　x 映射表 CV_32FC1/CV_32FC2
InputArray map2,　　　　　　y 映射表
int interpolation,　　　　　　　选择的插值方法，常见线性插值，可选择立方等
int borderMode,　　　　　　　BORDER_CONSTANT
const Scalar borderValue　　　color
)

头文件 quick_opencv.h：声明类与公共函数

#pragma once
#include <opencv2\opencv.hpp>
using namespace cv;class QuickDemo {
public:...void remap_Demo(Mat& image1);void MLS(Mat& src, std::vector<Point> p, std::vector<Point> q);void MLS(Mat& src, int* p, int* q, int rows, int cols);
};

主函数调用该类的公共成员函数

#include <opencv2\opencv.hpp>
#include <quick_opencv.h>
#include <iostream>
using namespace cv;int main(int argc, char** argv) {Mat src = imread("D:\\Desktop\\pandas_small22.png");if (src.empty()) {printf("Could not load images...\n");return -1;}QuickDemo qk;qk.remap_Demo(src);vector<Point> p{Point(30, 147), Point(147, 147), Point(268, 147), Point(112, 148),Point(186, 148), Point(98, 316), Point(211, 316)};vector<Point> q{ Point(28, 209), Point(126, 143), Point(282, 26), Point(71, 236), Point(136, 240), Point(79, 313), Point(190, 310)};qk.MLS(src1, p, q);int p_array[7][2] = { {30, 147}, {147, 147}, {268, 147}, {112, 148}, {186, 148}, {98, 316}, {211, 316} };int q_array[7][2] = { {28, 209}, {126, 143}, {282, 26},  {71, 236},  {136, 240}, {79, 313}, {190, 310} };qk.MLS(src1, (int *)p_array, (int*)q_array, 7, 2);waitKey(0);destroyAllWindows();return 0;
}

源文件 quick_demo.cpp：实现类与公共函数

void update_map(Mat& image, int index, Mat& x_map, Mat& y_map) {int height = image.rows;int width = image.cols;double h_41 = height * 0.25;double h_43 = height * 0.75;double w_41 = width * 0.25;double w_43 = width * 0.75;for (int h = 0; h < height; h++) {float* x_ptr = x_map.ptr<float>(h);float* y_ptr = y_map.ptr<float>(h);for (int w = 0; w < width; w++) {switch (index){case 0:if (h > h_41 && h < h_43 && w>w_41 && w < w_43) {*x_ptr++ = 2 * (w - w_41 + 0.5);*y_ptr++ = 2 * (h - h_41 + 0.5);}else{*x_ptr++ = 0;*y_ptr++ = 0;}break;case 1:*x_ptr++ = width - w - 1;*y_ptr++ = h;break;case 2:*x_ptr++ = w;*y_ptr++ = height - h - 1;break;case 3:*x_ptr++ = width - w - 1;*y_ptr++ = height - h - 1;break;}}}}void QuickDemo::remap_Demo(Mat& image) {Mat dst, x_map, y_map;int index = 0;x_map.create(image.size(), CV_32FC1);y_map.create(image.size(), CV_32FC1);int c = 0;while (true){c = waitKey(400);if ((char)c==27){break;}index = c % 4;update_map(image,index, x_map, y_map);remap(image, dst, x_map, y_map, INTER_LINEAR, BORDER_CONSTANT, Scalar(255, 0, 0));imshow("remap", dst);}
}

如上两个函数，update_map，用于更新remap的具体映射方法，remap_Demo为调用函数。

效果演示

运行程序，按键q退出；按键0，1，2，3，4 出现以下结果：

二、图像扭曲

MLS算法 图像扭曲 Image Deformation Using Moving Least Squares 论文。
最小二乘法（MLS）对图像进行变形 python 实现

整个算法公式：
$$f(v)=(v-p_{\ast }){\left(\sum _i{\widehat{p_i}}^T{w}_i\widehat{p_i}\right)}^{-1}\sum _j{\widehat{p_j}}^T{w}_j\widehat{p_j}+q_{\ast }\text{\hspace{1em}(1)}$$

其中: $v$ 是图像中的点，$p$ 是变形前的点，$q$ 是变形后的点，$w$ 是权重。

$$w_i=\frac{1}{|p_i-v{|}^{2\alpha }}\text{\hspace{1em}(2)}$$

$p_i$是图像中控制点的点集，$v$ 是原图像中的点，$\alpha$在程序中取1。

$$\begin{array}{rlrl}{p}_{\ast }=& \frac{{\sum }_i{w}_i{p}_i}{{\sum }_i{w}_i}（3）& \hat{p}=& p_i-p_{\ast }（5）\\ q_{\ast }=& \frac{{\sum }_i{w}_i{q}_i}{{\sum }_i{w}_i}（4）& \hat{q}=& q_i-q_{\ast }（6）\end{array}\text{\hspace{1em}(7)}$$

Point NewPoint(Point V, vector<Point> p, vector<Point> q){vector<float>W;Point p_star, q_star = Point(0, 0);for (int i = 0; i <= p.size() - 1; i++){float temp;if (p[i] == V){temp = INT_MAX;}else{temp = 1.0 / (((p[i].x - V.x) * (p[i].x - V.x)) + ((p[i].y - V.y) * (p[i].y - V.y)));}W.push_back(temp);}float px = 0, py = 0, qx = 0, qy = 0, W_sum = 0;for (int i = 0; i <= W.size() - 1; i++){px += W[i] * p[i].x;py += W[i] * p[i].y;qx += W[i] * q[i].x;qy += W[i] * q[i].y;W_sum += W[i];}p_star.x = px / W_sum;p_star.y = py / W_sum;q_star.x = qx / W_sum;q_star.y = qy / W_sum;vector<Point> p_hat, q_hat;for (int i = 0; i <= p.size() - 1; i++){p_hat.push_back(p[i] - p_star);q_hat.push_back(q[i] - q_star);}Mat pi_hat_t_ = Mat::zeros(2, 1, CV_32FC1);Mat_<float> pi_hat_t = pi_hat_t_;Mat pi_hat_ = Mat::zeros(1, 2, CV_32FC1);Mat_<float> pi_hat = pi_hat_;Mat M_1_ = Mat::zeros(2, 2, CV_32FC1);Mat_<float> M_1 = M_1_;for (int i = 0; i <= p_hat.size() - 1; i++){pi_hat_t.at<float>(0, 0) = p_hat[i].x;pi_hat_t.at<float>(1, 0) = p_hat[i].y;pi_hat.at<float>(0, 0) = p_hat[i].x;pi_hat.at<float>(0, 1) = p_hat[i].y;M_1 += pi_hat_t * W[i] * pi_hat;}Mat_<float> M_1_inv = M_1.inv();M_1 = M_1_inv;Mat pj_hat_t_ = Mat::zeros(2, 1, CV_32FC1);Mat_<float> pj_hat_t = pj_hat_t_;Mat qj_hat_ = Mat::zeros(1, 2, CV_32FC1);Mat_<float> qj_hat = qj_hat_;Mat M_2_ = Mat::zeros(2, 2, CV_32FC1);Mat_<float> M_2 = M_2_;for (int j = 0; j <= q.size() - 1; j++){pj_hat_t.at<float>(0, 0) = p_hat[j].x;pj_hat_t.at<float>(1, 0) = p_hat[j].y;qj_hat.at<float>(0, 0) = q_hat[j].x;qj_hat.at<float>(0, 1) = q_hat[j].y;M_2 += W[j] * pj_hat_t * qj_hat;}Mat_<float> M = M_1 * M_2;//ok//cout << "M = " << M << endl;Point x_p_star = V - p_star;Mat M_x_p_star_ = Mat::zeros(1, 2, CV_32FC1);Mat_<float> M_x_p_star = M_x_p_star_;M_x_p_star.at<float>(0, 0) = x_p_star.x;M_x_p_star.at<float>(0, 1) = x_p_star.y;Mat M_q_star_ = Mat::zeros(1, 2, CV_32FC1);Mat_<float> M_q_star = M_q_star_;M_q_star.at<float>(0, 0) = q_star.x;M_q_star.at<float>(0, 1) = q_star.y;Mat_<float> Lv = M_x_p_star * M + M_q_star;return Point(Lv.at<float>(0, 0), Lv.at<float>(0, 1));
}void QuickDemo::MLS(Mat& src, std::vector<Point> p, std::vector<Point> q){double time0 = static_cast<double>(getTickCount());Mat dst = Mat::zeros(src.rows, src.cols, CV_8UC3);for (int i = 0; i < src.rows; i++){for (int j = 0; j < src.cols; j++){Point old = Point(j, i);Point new_point = NewPoint(old, p, q);//cout << "old = " << old << "\tnew  = " << new_point << endl;dst.at<Vec3b>(i, j) = src.at<Vec3b>(abs(new_point.y), abs(new_point.x));}}double time1 = static_cast<double>(getTickCount());cout << "Total cost time is " << ((time1 - time0) / getTickFrequency()) << "seconds" << endl;imshow("dst_msl", dst);
}

重载函数

Point NewPoint(Point V, float* W, int* p, int* q , float* p_hat, float* q_hat, int rows, int cols) {Point p_star, q_star = Point(0, 0);float temp = 0;float px = 0, py = 0, qx = 0, qy = 0, W_sum = 0;for (int i = 0; i < rows; i++) {int p_0 = *(p + i * cols);int p_1 = *(p + i * cols + 1);if (!(p_0 == V.x && p_1 == V.y)) {temp = 1.0 / (((p_0 - V.x) * (p_0 - V.x)) + ((p_1 - V.y) * (p_1 - V.y)));}else {temp = INT_MAX;}W[i] = temp;px += temp * p_0;py += temp * p_1;qx += temp * (*(q + i * cols));qy += temp * (*(q + i * cols + 1));W_sum += temp;}p_star.x = px / W_sum;p_star.y = py / W_sum;q_star.x = qx / W_sum;q_star.y = qy / W_sum;for (int i = 0; i < rows; i++) {*(p_hat + i * cols) = *(p + i * cols) - p_star.x;*(p_hat + i * cols + 1) = *(p + i * cols + 1) - p_star.y;*(q_hat + i * cols) = *(q + i * cols) - p_star.x;*(q_hat + i * cols + 1) = *(q + i * cols + 1) - p_star.y;}// ====================================Mat pi_hat_t_ = Mat::zeros(2, 1, CV_32FC1);Mat_<float> pi_hat_t = pi_hat_t_;Mat pi_hat_ = Mat::zeros(1, 2, CV_32FC1);Mat_<float> pi_hat = pi_hat_;Mat M_1_ = Mat::zeros(2, 2, CV_32FC1);Mat_<float> M_1 = M_1_;// ====================================Mat pj_hat_t_ = Mat::zeros(2, 1, CV_32FC1);Mat_<float> pj_hat_t = pj_hat_t_;Mat qj_hat_ = Mat::zeros(1, 2, CV_32FC1);Mat_<float> qj_hat = qj_hat_;Mat M_2_ = Mat::zeros(2, 2, CV_32FC1);Mat_<float> M_2 = M_2_;// ====================================for (int i = 0; i < rows; i++) {float p_hat_x = *(p_hat + i * cols);float p_hat_y = *(p_hat + i * cols + 1);pi_hat_t.at<float>(0, 0) = p_hat_x;pi_hat_t.at<float>(1, 0) = p_hat_y;pi_hat.at<float>(0, 0) = p_hat_x;pi_hat.at<float>(0, 1) = p_hat_y;M_1 += pi_hat_t * W[i] * pi_hat;pj_hat_t.at<float>(0, 0) = p_hat_x;pj_hat_t.at<float>(1, 0) = p_hat_y;qj_hat.at<float>(0, 0) = *(q_hat + i * cols);qj_hat.at<float>(0, 1) = *(q_hat + i * cols + 1);M_2 += pj_hat_t * W[i] * qj_hat;}Mat_<float> M_1_inv = M_1.inv();M_1 = M_1_inv;Mat_<float> M = M_1 * M_2;//=====================================//// 	  如下为总公式计算////======================================Point x_p_star = V - p_star;Mat M_x_p_star_ = Mat::zeros(1, 2, CV_32FC1);Mat_<float> M_x_p_star = M_x_p_star_;M_x_p_star.at<float>(0, 0) = x_p_star.x;M_x_p_star.at<float>(0, 1) = x_p_star.y;Mat M_q_star_ = Mat::zeros(1, 2, CV_32FC1);Mat_<float> M_q_star = M_q_star_;M_q_star.at<float>(0, 0) = q_star.x;M_q_star.at<float>(0, 1) = q_star.y;Mat_<float> Lv = M_x_p_star * M + M_q_star;return Point(Lv.at<float>(0, 0), Lv.at<float>(0, 1));}void QuickDemo::MLS(Mat& src, int* p, int* q, int rows, int cols) {double time0 = static_cast<double>(getTickCount());Mat dst = Mat::zeros(src.rows, src.cols, CV_8UC3);assert(7 == rows);               // 若断言失败请修改如下三个数组的长度为rowsfloat W[7] = { 0 };              // 权重长度为p数组长度：rows=7float p_hat[7][2] = { 0 };       // p_hat长度为p数组长度：rows=7float q_hat[7][2] = { 0 };       // q_hat长度为p数组长度：rows=7for (int i = 0; i < src.rows; i++) {for (int j = 0; j < src.cols; j++) {Point new_point = NewPoint(Point(j, i), W, p, q, (float*)p_hat, (float*)p_hat, rows, cols);//cout << "old = " << old << "\tnew  = " << new_point << endl;dst.at<Vec3b>(i, j) = src.at<Vec3b>(abs(new_point.y), abs(new_point.x));//cout << "src.at<uchar> = " << src.at<Vec3b>(new_point.y,new_point.x) << endl;}}double time1 = static_cast<double>(getTickCount());cout << "Total cost time is " << ((time1 - time0) / getTickFrequency()) << "seconds" << endl;imshow("dst_msl", dst);
}

重载后，就快了100~180ms，然并卵。
效果图：

鸣谢与拓展阅读：
使用范例 记录四图像处理之瘦脸 MLS算法 C++实现
OpenCV局部变形算法探究
基于移动最小二乘（MLS）的图像扭曲刚性变形python实现
使用重映射实现图像的局部扭曲 来实现 图像增强。