C# 实现TPS薄板样条差值

先是实现了一个类

using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;
using Accord;

namespace QRDetectAndDecode
{
    class ThinPlateSpline
    {
        private int n;
        private double[,] K;
        private double[,] P;
        private double[,] L;
        private double[] Vx, Vy;
        private double[] paramsX, paramsY;

        public ThinPlateSpline(Vector2[] src, Vector2[] dst)
        {
            n = src.Length;
            K = new double[n, n];
            P = new double[n, 3];
            L = new double[n + 3, n + 3];
            Vx = new double[n + 3];
            Vy = new double[n + 3];

            // 计算K
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    if (i == j)
                        K[i, j] = 0;
                    else
                    {
                        double r2 = Math.Pow(src[i].X - src[j].X, 2) + Math.Pow(src[i].Y - src[j].Y, 2);
                        K[i, j] = r2 == 0 ? 0 : r2 * Math.Log(r2 + 1e-20);
                    }
                }
            }
            // 计算P
            for (int i = 0; i < n; i++)
            {
                P[i, 0] = 1;
                P[i, 1] = src[i].X;
                P[i, 2] = src[i].Y;
            }
            // 组装L
            for (int i = 0; i < n; i++)
                for (int j = 0; j < n; j++)
                    L[i, j] = K[i, j];
            for (int i = 0; i < n; i++)
                for (int j = 0; j < 3; j++)
                {
                    L[i, n + j] = P[i, j];
                    L[n + j, i] = P[i, j];
                }
            // 填充右下角0
            for (int i = n; i < n + 3; i++)
                for (int j = n; j < n + 3; j++)
                    L[i, j] = 0;

            // 组装V
            for (int i = 0; i < n; i++)
            {
                Vx[i] = dst[i].X;
                Vy[i] = dst[i].Y;
            }
            for (int i = n; i < n + 3; i++)
            {
                Vx[i] = 0;
                Vy[i] = 0;
            }

            // 解线性方程组
            paramsX = SolveLinearSystem(L, Vx);
            paramsY = SolveLinearSystem(L, Vy);
        }

        // TPS映射单点
        public Vector2 Transform(Vector2 pt)
        {
            double[] U = new double[n];
            for (int i = 0; i < n; i++)
            {
                double r2 = Math.Pow(pt.X - srcPts[i].X, 2) + Math.Pow(pt.Y - srcPts[i].Y, 2);
                U[i] = r2 == 0 ? 0 : r2 * Math.Log(r2 + 1e-20);
            }
            double[] P = new double[] { 1, pt.X, pt.Y };

            double x = 0, y = 0;
            for (int i = 0; i < n; i++)
            {
                x += paramsX[i] * U[i];
                y += paramsY[i] * U[i];
            }
            for (int i = 0; i < 3; i++)
            {
                x += paramsX[n + i] * P[i];
                y += paramsY[n + i] * P[i];
            }
            return new Vector2((float)x, (float)y);
        }

        // 解线性方程组Ax = b
        private double[] SolveLinearSystem(double[,] A, double[] b)
        {
            int n = b.Length;
            var mat = Accord.Math.Matrix.MemberwiseClone(A); // 需要 Accord.Math
            var vec = b.ToArray();
            return Accord.Math.Matrix.Solve(mat, vec, leastSquares: false);
        }

        private Vector2[] srcPts;
        public void SetSourcePoints(Vector2[] src)
        {
            srcPts = src;
        }
    }
}

下面是这个类的调用方法

        public static Mat PPTPS(Point2f[] pts, Mat srcImg,int size)
        {
            Vector2[] src6 = new Vector2[6];
            for (int i = 0; i < 4; i++)
                src6[i] = new Vector2(pts[i].X, pts[i].Y);

            // 2. 计算上下中点
            src6[4] = new Vector2((pts[0].X + pts[1].X) / 2, (pts[0].Y + pts[1].Y) / 2); // 上中点
            src6[5] = new Vector2((pts[2].X + pts[3].X) / 2, (pts[2].Y + pts[3].Y) / 2); // 下中点
            src6[4].Y += 2;
            src6[5].Y += 2;

            Vector2[] dst6 = new Vector2[6];

            var XX = src6[0].X;
            var YY = src6[0].Y;
            dst6[0] = new Vector2(XX, YY);                   // 左上
            dst6[1] = new Vector2(XX + size*2, YY);            // 右上
            dst6[2] = new Vector2(XX + size * 2, YY + size * 2);    // 右下
            dst6[3] = new Vector2(XX, YY + size * 2);           // 左下
            dst6[4] = new Vector2(XX + size, YY);        // 上中点
            dst6[5] = new Vector2(XX + size, YY + size * 2);// 下中点

            var tps = new ThinPlateSpline(dst6, src6);
            tps.SetSourcePoints(src6);

            Mat tpsResult = new(srcImg.Height, srcImg.Width, srcImg.Type());

            for (int y = 0; y < srcImg.Height; y++)
            {
                for (int x = 0; x < srcImg.Width; x++)
                {
                    Vector2 dstPt = new(x, y);
                    // 反向映射到原图
                    Vector2 srcPt = tps.Transform(dstPt);

                    // 边界检查
                    int srcX = (int)Math.Round(srcPt.X);
                    int srcY = (int)Math.Round(srcPt.Y);
                    if (srcX >= 0 && srcX < srcImg.Cols && srcY >= 0 && srcY < srcImg.Rows)
                    {
                        tpsResult.Set(y, x, srcImg.At<Vec3b>(srcY, srcX));
                    }
                    else
                    {
                        tpsResult.Set(y, x, new Vec3b(0, 0, 0)); // 超出边界设为黑色
                    }
                }
            }
            return tpsResult;
        }

            src6[4].Y += 2;
            src6[5].Y += 2;

这两句是我用来调整图像形变效果用的

后续还要优化TPS的插值方法，现在图像效果不够好