TPS进行轮廓变形

#include <opencv2/opencv.hpp>
//#include <opencv2/xfeatures2d.hpp>
#include <opencv2/features2d.hpp>
#include <vector>

using namespace cv;
using namespace std;
//using namespace cv::xfeatures2d;

class ThinPlateSpline {
public:
	void compute(const vector<Point2f>& source, const vector<Point2f>& target) {
		int n = source.size();
		if (n < 3) return;

		Mat L = Mat::zeros(n + 3, n + 3, CV_32F);
		Mat Vx = Mat::zeros(n + 3, 1, CV_32F);
		Mat Vy = Mat::zeros(n + 3, 1, CV_32F);

		// 构建K矩阵
		for (int i = 0; i < n; ++i) {
			for (int j = 0; j < n; ++j) {
				float r = norm(source[i] - source[j]);
				L.at<float>(i, j) = r * r * log(r + 1e-6);
			}
		}

		// 填充P矩阵及其转置
		for (int i = 0; i < n; ++i) {
			L.at<float>(i, n) = 1.0;
			L.at<float>(i, n + 1) = source[i].x;
			L.at<float>(i, n + 2) = source[i].y;

			L.at<float>(n, i) = 1.0;
			L.at<float>(n + 1, i) = source[i].x;
			L.at<float>(n + 2, i) = source[i].y;
		}

		// 设置V向量为目标坐标
		for (int i = 0; i < n; ++i) {
			Vx.at<float>(i) = target[i].x;
			Vy.at<float>(i) = target[i].y;
		}

		// 解方程获取系数
		solve(L, Vx, Wx, DECOMP_SVD);
		solve(L, Vy, Wy, DECOMP_SVD); //DECOMP_LU

		src_points = source;
	}

	Point2f warpPoint(const Point2f& p) const {
		float fx = Wx.at<float>(src_points.size()) + Wx.at<float>(src_points.size() + 1) * p.x + Wx.at<float>(src_points.size() + 2) * p.y;
		float fy = Wy.at<float>(src_points.size()) + Wy.at<float>(src_points.size() + 1) * p.x + Wy.at<float>(src_points.size() + 2) * p.y;

		for (int i = 0; i < src_points.size(); ++i) {
			float r = norm(p - src_points[i]);
			float U = r * r * log(r + 1e-6);
			fx += Wx.at<float>(i) * U;
			fy += Wy.at<float>(i) * U;
		}

		return Point2f(fx, fy);
	}

private:
	Mat Wx, Wy;
	vector<Point2f> src_points;
};

int main() {
	// 读取图像
	Mat img_template = imread("a.jpg", IMREAD_GRAYSCALE);
	Mat img_actual = imread("b.jpg", IMREAD_GRAYSCALE);
	//Mat img_template = imread("1.png", IMREAD_GRAYSCALE);
	//Mat img_actual = imread("2.png", IMREAD_GRAYSCALE);
	cv::resize(img_actual, img_actual, img_template.size());

	// SIFT特征检测
	Ptr<SIFT> sift = SIFT::create();
	vector<KeyPoint> kp_template, kp_actual;
	Mat desc_template, desc_actual;
	sift->detectAndCompute(img_template, noArray(), kp_template, desc_template);
	sift->detectAndCompute(img_actual, noArray(), kp_actual, desc_actual);

	// KNN匹配
	BFMatcher matcher(NORM_L2);
	vector<vector<DMatch>> knn_matches;
	matcher.knnMatch(desc_template, desc_actual, knn_matches, 2);

	// 筛选匹配
	vector<DMatch> good_matches;
	for (size_t i = 0; i < knn_matches.size(); i++) {
		if (knn_matches[i][0].distance < 0.7 * knn_matches[i][1].distance) {
			good_matches.push_back(knn_matches[i][0]);
		}
	}

	// 获取对应点
	vector<Point2f> src_pts, dst_pts;
	for (auto& m : good_matches) {
		src_pts.push_back(kp_template[m.queryIdx].pt);
		dst_pts.push_back(kp_actual[m.trainIdx].pt);
	}

	// RANSAC筛选
	Mat mask;
	findHomography(src_pts, dst_pts, RANSAC, 3, mask);
	vector<Point2f> reliable_src, reliable_dst;
	for (int i = 0; i < mask.rows; ++i) {
		if (mask.at<uchar>(i)) {
			reliable_src.push_back(src_pts[i]);
			reliable_dst.push_back(dst_pts[i]);
		}
	}

	if (reliable_src.size() < 3) {
		cerr << "Not enough points." << endl;
		return -1;
	}

	// 计算TPS
	ThinPlateSpline tps;
	tps.compute(reliable_src, reliable_dst);

	// 提取模板轮廓
	cv:Mat binImage;
	vector<vector<Point>> contours;
	cv::threshold(img_template, binImage, 50, 255, cv::THRESH_BINARY); //THRESH_BINARY_INV
	findContours(binImage, contours, RETR_EXTERNAL, CHAIN_APPROX_SIMPLE);

	// 应用TPS变形
	vector<vector<Point2f>> warped_contours;
	for (auto& contour : contours) {
		vector<Point2f> warped_contour;
		for (auto& pt : contour) {
			Point2f warped_pt = tps.warpPoint(Point2f(pt));
			warped_contour.push_back(warped_pt);
		}
		warped_contours.push_back(warped_contour);
	}

	// 绘制结果
	Mat img_color;
	cvtColor(img_actual, img_color, COLOR_GRAY2BGR);
	drawContours(img_color, contours, -1, Scalar(0, 255, 0), 1);

	for (size_t m = 0; m < src_pts.size(); m++) {
		cv::line(img_color, src_pts[m], dst_pts[m], Scalar(0, 0, 255), 1);
	}

	for (auto& contour : warped_contours) {
		vector<Point> int_contour;
		for (auto& pt : contour) {
			int_contour.push_back(Point(cvRound(pt.x), cvRound(pt.y)));
		}
		drawContours(img_color, vector<vector<Point>>{int_contour}, -1, Scalar(255, 0, 0), 1);
	}

	imshow("Result", img_color);
	waitKey();

	return 0;
}