[数据结构与算法]C#，贝塞耳插值（Bessel‘s interpolation）的算法与源代码

Friedrich Wilhelm Bessel

首先要区别于另外一个读音接近的插值算法：贝塞尔插值（Bézier）。

（1）读音接近，但不是一个人；

（2）一个是多项式（整体）插值，一个是分段插值；

（3）一个已经很少用，一个还是应用主力；

贝塞耳插值（Bessel's interpolation）是一种等距节点插值方法，适用于被插值节点z位于插值区间中部且位于两相邻插值点的中点附近的情况。

using System;
using System.Text;
using System.Collections;
using System.Collections.Generic;

namespace Legalsoft.Truffer.Algorithm
{
    public partial class TPoint
    {
        public double X { get; set; } = 0.0;
        public double Y { get; set; } = 0.0;
        public double Z { get; set; } = 0.0;
        public TPoint()
        {
        }
        public TPoint(double x, double y)
        {
            X = x; Y = y;
        }
        public TPoint(double x, double y, double z)
        {
            X = x; Y = y; Z = z;
        }

        public double Distance(TPoint p1)
        {
            double ds = (p1.X - this.X) * (p1.X - this.X) + (p1.Y - this.Y) * (p1.Y - this.Y);
            if (ds <= float.Epsilon) return 0.0;
            return Math.Sqrt(ds);
        }

        public static double Distance(TPoint p1, TPoint p2)
        {
            double ds = (p1.X - p2.X) * (p1.X - p2.X) + (p1.Y - p2.Y) * (p1.Y - p2.Y);
            if (ds <= float.Epsilon) return 0.0;
            return Math.Sqrt(ds);
        }
    }

    public static partial class Algorithm_Gallery
    {
        private static double U_Calculate(double u, int n)
        {
            if (n == 0)
            {
                return 1.0;
            }
            double temp = u;
            for (int i = 1; i <= n / 2; i++)
            {
                temp = temp * (u - i);
            }
            for (int i = 1; i < n / 2; i++)
            {
                temp = temp * (u + i);
            }
            return temp;
        }

        private static int Fact(int n)
        {
            int f = 1;
            for (int i = 2; i <= n; i++)
            {
                f *= i;
            }
            return f;
        }

        public static double Bessel_Interpolation(List<TPoint> points, double value)
        {
            int n = points.Count;
            double[,] y = new double[n, n];
            for (int i = 0; i < n; i++)
            {
                y[i, 0] = points[i].Y;
            }
            for (int i = 1; i < n; i++)
            {
                for (int j = 0; j < n - i; j++)
                {
                    y[j, i] = y[j + 1, i - 1] - y[j, i - 1];
                }
            }

            double sum = (y[2, 0] + y[3, 0]) / 2;

            int k;
            if ((n % 2) > 0)
            {
                k = n / 2;
            }
            else
            {
                k = (n / 2) - 1; // origin for even
            }
            double u = (value - points[k].X) / (points[1].X - points[0].X);

            for (int i = 1; i < n; i++)
            {
                if ((i % 2) > 0)
                {
                    sum = sum + ((u - 0.5) * U_Calculate(u, i - 1) * y[k, i]) / Fact(i);
                }
                else
                {
                    sum = sum + (U_Calculate(u, i) * (y[k, i] + y[--k, i]) / (Fact(i) * 2));
                }
            }
            return sum;
        }
    }
}

创作打卡挑战赛

赢取流量/现金/CSDN周边激励大奖