3次スプライン補間基本クラス - わかばテクノロジ

補間
3次スプライン補間基本クラス; Cubic Spline

概要

スプライン補間とは\(N\)個の点の間の区間を多項式を用いて補間する手法のことで、補間する曲線をスプライン曲線という。
3次スプライン補間は以下の3次多項式を用いる。

\(\quad f(x)=a+b x + c^2 x + d x^3\)

この\(a, b, c, d\)は区間の左端\(x = 0\)における制御点の値と傾き\(v_0, g_0\)と区間の右端\(x = 1\)における制御点の値と傾き\(v_1, g_1\)によって一意に定まる。

\(\quad f(0) = v_0 = a, \ f'(0) = g_0 = b, \ f(1) = v_1 = a + b + c + d, \ f'(1) = g_1 = b + 2 c + 3 d \)
\(\quad \rightarrow a = v_0, \ b = g_0, \ c = -3 (v_0 - v1) - 2 g_0 - g_1, \ d = 2(v_0 - v_1) + g_0 + g_1 \)

制御点における傾きの算出方法はスプライン補間の手法によって様々あり、用途によって選択する必要がある。
当サイトでは現在、秋間スプライン、単調スプライン、Catmull-Romスプラインを掲載している。
この記事ではそれらの共通する機能となる基本クラスを扱う。

スプライン補間構成

・スプライン基本クラス

　　・3次スプライン基本クラス

　　　　・制御点とその近傍2点から傾きを算出する3次スプライン基本クラス

　　　　　　・単調スプライン

　　　　　　・Catmull-Romスプライン

　　　　・制御点とその近傍4点から傾きを算出する3次スプライン基本クラス

　　　　　　・秋間スプライン

・終端タイプ列挙型

・多次元スプライン補間ジェネリッククラス

スプライン基本クラスソースコード


namespace SplineInterpolation {

    /// <summary>スプライン基本クラス</summary>
    public abstract class Spline {

        /// <summary>制御点</summary>
        protected List<double> v;

        /// <summary>コンストラクタ</summary>
        public Spline(EndType type = EndType.Open) {
            this.EndType = type;
            this.v = new List<double>();
        }

        /// <summary>終端タイプ</summary>
        public EndType EndType { get; private set; }

        /// <summary>制御点数</summary>
        public int Points => v.Count;

        /// <summary>制御点を設定</summary>
        public void Set(params double[] v) {
            Set(v, v.Length);
        }

        /// <summary>制御点を設定</summary>
        public abstract void Set(double[] v, int points);

        /// <summary>制御点を挿入</summary>
        public abstract void Insert(int index, double new_v);

        /// <summary>制御点を除去</summary>
        public abstract void Remove(int index);

        /// <summary>初期化</summary>
        public virtual void Initialize() {
            v.Clear();
        }

        /// <summary>補間値</summary>
        public abstract double Value(double t);

        /// <summary>補間微分値</summary>
        public abstract double Diff(double t);

        /// <summary>補間N次微分値</summary>
        public abstract double Diff(double t, uint n);

        /// <summary>制御点番号が有効であるか判定</summary>
        public bool IsInRange(int index) {
            return (index >= 0) || (index < Points);
        }

        /// <summary>制御点の値を取得</summary>
        public double GetPoint(int index) {
            return v[index];
        }

        /// <summary>制御点の値を設定</summary>
        public abstract void SetPoint(int index, double set_v);

        /// <summary>等しいか判定</summary>
        public static bool operator ==(Spline s1, Spline s2) {
            if(ReferenceEquals(s1, s2)) {
                return true;
            }
            if((object)s1 == null || (object)s2 == null) {
                return false;
            }
            if(s1.EndType != s2.EndType) {
                return false;
            }

            return s1.v.SequenceEqual(s2.v);
        }

        /// <summary>異なるか判定</summary>
        public static bool operator !=(Spline s1, Spline s2) {
            return !(s1 == s2);
        }

        /// <summary>等しいか判定</summary>
        public override bool Equals(object obj) {
            return obj is Spline ? (Spline)obj == this : false;
        }

        /// <summary>ハッシュ値</summary>
        public override int GetHashCode() {
            return Points > 0 ? v[0].GetHashCode() : 0;
        }

        /// <summary>コンストラクタ ジェネリッククラス用</summary>
        public static SplineType Construct<SplineType>(EndType type = EndType.Open) {
            return (SplineType)typeof(SplineType).GetConstructor(new Type[] { typeof(EndType) }).Invoke(new object[] { type });
        }
    }
}

3次スプライン基本クラスソースコード


namespace SplineInterpolation {

    /// <summary>3次スプライン補間基本クラス</summary>
    public abstract class CubicSpline : Spline {

        /// <summary>区間リスト</summary>
        protected List<CubicSegment> segment_list;

        /// <summary>コンストラクタ</summary>
        public CubicSpline(EndType type = EndType.Open) : base(type) {
            this.segment_list = new List<CubicSegment>();
        }

        /// <summary>制御点を設定</summary>
        public sealed override void Set(double[] v, int points) {
            Initialize();

            if(points == 0) {
                return;
            }

            if(v == null || v.Length <= 0 || v.Length < points) {
                throw new ArgumentException();
            }

            this.v.AddRange(v.Take(points));

            segment_list = new List<CubicSegment>(points);

            for(int i = 0; i < points; i++) {
                segment_list.Add(new CubicSegment());
            }

            for(int i = 0; i < Points; i++) {
                ReflashSegment(i);
            }
        }

        /// <summary>初期化</summary>
        public sealed override void Initialize() {
            base.Initialize();
            segment_list.Clear();
        }

        /// <summary>補間値</summary>
        public sealed override double Value(double t) {
            if(Points <= 0 || double.IsNaN(t) || double.IsInfinity(t)) {
                return double.NaN;
            }

            if(EndType == EndType.Open) {
                int pos = t < (Points - 1) ? (int)t : (Points - 1);
                if(t >= 0) {
                    double h = t - pos;
                    return segment_list[pos].Value(h);
                }
                return segment_list[0].a + t * segment_list[0].b;
            }
            else {
                double h = t - Math.Floor(t);
                int pos = ((int)(Math.Floor(t)) % Points + Points) % Points;
                return segment_list[pos].Value(h);
            }
        }

        /// <summary>補間微分値</summary>
        public sealed override double Diff(double t) {
            if(Points <= 0 || double.IsNaN(t) || double.IsInfinity(t)) {
                return double.NaN;
            }

            if(EndType == EndType.Open) {
                int pos = t < (Points - 1) ? (int)t : (Points - 1);
                if(t >= 0) {
                    double h = t - pos;
                    return segment_list[pos].Diff(h);
                }
                return segment_list[0].b;
            }
            else {
                double h = t - Math.Floor(t);
                int pos = ((int)(Math.Floor(t)) % Points + Points) % Points;
                return segment_list[pos].Diff(h);
            }
        }

        /// <summary>補間N次微分値</summary>
        public sealed override double Diff(double t, uint n) {
            if(Points <= 0 || double.IsNaN(t) || double.IsInfinity(t)) {
                return double.NaN;
            }

            if(EndType == EndType.Open) {
                if(n == 0)
                    return Value(t);
                if(n == 1)
                    return Diff(t);
                if(n > 3)
                    return 0;

                int pos = t < (Points - 1) ? (int)t : (Points - 1);
                if(t >= 0) {
                    double h = t - pos;
                    return segment_list[pos].Diff(h, n);
                }
                return 0;
            }
            else {
                double h = t - Math.Floor(t);
                int pos = ((int)(Math.Floor(t)) % Points + Points) % Points;
                return segment_list[pos].Diff(h, n);
            }
        }

        /// <summary>制御点における傾き</summary>
        protected abstract double Grad(int index);

        /// <summary>区間更新</summary>
        protected void ReflashSegment(int index) {
            if(Points <= 1) {
                segment_list[0] = CubicSegment.Constant(v[0]);
                return;
            }

            if(EndType == EndType.Close) {
                index = (index % Points + Points) % Points;
            }

            if(index < 0 || index >= Points) {
                return;
            }

            if(EndType == EndType.Open) {
                if(index == 0) {
                    segment_list[0] = CubicSegment.LeftEnd(v[0], v[1], Grad(1));
                }
                else if(index >= Points - 2) {
                    segment_list[Points - 2] = CubicSegment.RightEnd(v[Points - 2], v[Points - 1], Grad(Points - 2));

                    double g = segment_list[Points - 2].b + 2 * segment_list[Points - 2].c;
                    segment_list[Points - 1] = CubicSegment.Linear(v[Points - 1], v[Points - 1] + g);
                }
                else {
                    segment_list[index] = CubicSegment.Interval(v[index], v[index + 1], Grad(index), Grad(index + 1));
                }
            }
            else {
                if(index == Points - 1) {
                    segment_list[Points - 1] = CubicSegment.Interval(v[Points - 1], v[0], Grad(Points - 1), Grad(0));
                }
                else {
                    segment_list[index] = CubicSegment.Interval(v[index], v[index + 1], Grad(index), Grad(index + 1));
                }
            }
        }

        /// <summary>等しいか判定</summary>
        public static bool operator ==(CubicSpline s1, CubicSpline s2) {
            if((Spline)s1 != (Spline)s2) {
                return false;
            }

            for(int i = 0, Points = s1.segment_list.Count; i < Points; i++) {
                if(s1.segment_list[i] != s2.segment_list[i]) {
                    return false;
                }
            }
            return true;
        }

        /// <summary>異なるか判定</summary>
        public static bool operator !=(CubicSpline s1, CubicSpline s2) {
            return !(s1 == s2);
        }

        /// <summary>等しいか判定</summary>
        public override bool Equals(object obj) {
            return obj is CubicSpline ? (CubicSpline)obj == this : false;
        }

        /// <summary>ハッシュ値</summary>
        public override int GetHashCode() {
            return Points > 0 ? segment_list[0].GetHashCode() : 0;
        }
    }
}

制御点とその近傍2点から傾きを算出する3次スプライン基本クラスソースコード


namespace SplineInterpolation {

    /// <summary>制御点とその近傍2点から傾きを算出する3次スプライン基本クラス</summary>
    public abstract class CubicSplineNeighbor2 : CubicSpline {

        /// <summary>コンストラクタ</summary>
        public CubicSplineNeighbor2(EndType type = EndType.Open) : base(type) {
        }

        /// <summary>制御点を挿入</summary>
        public sealed override void Insert(int index, double new_v) {
            if(index < 0 || index > Points) {
                throw new ArgumentException();
            }

            v.Insert(index, new_v);
            segment_list.Insert(index, new CubicSegment());

            if(Points > 4) {
                for(int i = index - 2; i <= index + 1; i++) {
                    ReflashSegment(i);
                }
            }
            else {
                for(int i = 0; i < Points; i++) {
                    ReflashSegment(i);
                }
            }
        }

        /// <summary>制御点を除去</summary>
        public sealed override void Remove(int index) {
            if(index < 0 || index >= Points) {
                throw new ArgumentException();
            }

            v.RemoveAt(index);
            segment_list.RemoveAt(index);

            if(Points > 3) {
                for(int i = index - 2; i <= index; i++) {
                    ReflashSegment(i);
                }
            }
            else {
                for(int i = 0; i < Points; i++) {
                    ReflashSegment(i);
                }
            }
        }

        /// <summary>制御点を設定</summary>
        public sealed override void SetPoint(int index, double set_v) {
            v[index] = set_v;

            if(Points > 4) {
                for(int i = index - 2; i <= index + 1; i++) {
                    ReflashSegment(i);
                }
            }
            else {
                for(int i = 0; i < Points; i++) {
                    ReflashSegment(i);
                }
            }
        }

        /// <summary>制御点における傾き</summary>
        protected abstract double Grad(double vm1, double v0, double vp1);

        /// <summary>制御点における傾き</summary>
        protected sealed override double Grad(int index) {
            if(index < 0 || index >= Points) {
                throw new ArgumentException();
            }

            if(Points <= 1) {
                return 0;
            }

            if(EndType == EndType.Open) {
                if(Points <= 2) {
                    return v[1] - v[0];
                }

                if(index == 0) {
                    return Grad(2 * v[0] - v[1], v[0], v[1]);
                }
                if(index == Points - 1) {
                    return Grad(v[Points - 2], v[Points - 1], 2 * v[Points - 1] - v[Points - 2]);
                }
            }
            else {
                if(Points <= 2) {
                    return 0;
                }

                if(index == 0) {
                    return Grad(v[Points - 1], v[0], v[1]);
                }
                if(index == Points - 1) {
                    return Grad(v[Points - 2], v[Points - 1], v[0]);
                }
            }
            return Grad(v[index - 1], v[index], v[index + 1]);
        }
    }
}

制御点とその近傍4点から傾きを算出する3次スプライン基本クラスソースコード


namespace SplineInterpolation {

    /// <summary>制御点とその近傍4点から傾きを算出する3次スプライン基本クラス</summary>
    public abstract class CubicSplineNeighbor4 : CubicSpline {

        /// <summary>コンストラクタ</summary>
        public CubicSplineNeighbor4(EndType type = EndType.Open) : base(type) {
        }

        /// <summary>制御点を挿入</summary>
        public sealed override void Insert(int index, double new_v) {
            if(index < 0 || index > Points) {
                throw new ArgumentException();
            }

            v.Insert(index, new_v);
            segment_list.Insert(index, new CubicSegment());

            if(Points > 6) {
                for(int i = index - 3; i <= index + 2; i++) {
                    ReflashSegment(i);
                }
            }
            else {
                for(int i = 0; i < Points; i++) {
                    ReflashSegment(i);
                }
            }
        }

        /// <summary>制御点を除去</summary>
        public sealed override void Remove(int index) {
            if(index < 0 || index >= Points) {
                throw new ArgumentException();
            }

            v.RemoveAt(index);
            segment_list.RemoveAt(index);

            if(Points > 5) {
                for(int i = index - 3; i <= index + 1; i++) {
                    ReflashSegment(i);
                }
            }
            else {
                for(int i = 0; i < Points; i++) {
                    ReflashSegment(i);
                }
            }
        }

        /// <summary>制御点を設定</summary>
        public sealed override void SetPoint(int index, double set_v) {
            v[index] = set_v;

            if(Points > 6) {
                for(int i = index - 3; i <= index + 2; i++) {
                    ReflashSegment(i);
                }
            }
            else {
                for(int i = 0; i < Points; i++) {
                    ReflashSegment(i);
                }
            }
        }

        /// <summary>制御点における傾き</summary>
        protected abstract double Grad(double vm2, double vm1, double v0, double vp1, double vp2);

        /// <summary>制御点における傾き</summary>
        protected sealed override double Grad(int index) {
            if(index < 0 || index >= Points) {
                throw new ArgumentException();
            }

            if(Points <= 1) {
                return 0;
            }

            if(EndType == EndType.Open) {
                if(Points <= 2) {
                    return v[1] - v[0];
                }

                if(Points <= 3 && index == 1) {
                    return Grad(2 * v[0] - v[1], v[0], v[1], v[2], 2 * v[2] - v[1]);
                }

                if(index == 0) {
                    return Grad(3 * v[0] - 2 * v[1], 2 * v[0] - v[1], v[0], v[1], v[2]);
                }
                if(index == 1) {
                    return Grad(2 * v[0] - v[1], v[0], v[1], v[2], v[3]);
                }
                if(index == Points - 2) {
                    return Grad(v[Points - 4], v[Points - 3], v[Points - 2], v[Points - 1], 2 * v[Points - 1] - v[Points - 2]);
                }
                if(index == Points - 1) {
                    return Grad(v[Points - 3], v[Points - 2], v[Points - 1], 2 * v[Points - 1] - v[Points - 2], 3 * v[Points - 1] - 2 * v[Points - 2]);
                }
            }
            else {
                if(Points <= 2) {
                    return 0;
                }

                if(Points <= 3) {
                    switch(index) {
                        case 0:
                            return Grad(v[1], v[2], v[0], v[1], v[2]);
                        case 1:
                            return Grad(v[2], v[0], v[1], v[2], v[0]);
                        default:
                            return Grad(v[0], v[1], v[2], v[0], v[1]);
                    }
                }

                if(index == 0) {
                    return Grad(v[Points - 2], v[Points - 1], v[0], v[1], v[2]);
                }
                if(index == 1) {
                    return Grad(v[Points - 1], v[0], v[1], v[2], v[3]);
                }
                if(index == Points - 1) {
                    return Grad(v[Points - 3], v[Points - 2], v[Points - 1], v[0], v[1]);
                }
                if(index == Points - 2) {
                    return Grad(v[Points - 4], v[Points - 3], v[Points - 2], v[Points - 1], v[0]);
                }

            }
            return Grad(v[index - 2], v[index - 1], v[index], v[index + 1], v[index + 2]);
        }
    }
}

終端タイプ列挙型


namespace SplineInterpolation {

    /// <summary>終端タイプ</summary>
    public enum EndType {
        /// <summary>終端をループさせない</summary>
        Open,
        /// <summary>終端をループさせる</summary>
        Close
    };
}

多次元スプライン補間ジェネリッククラス


namespace SplineInterpolation {

    /// <summary>多次元スプライン補間ジェネリッククラス</summary>
    public class SplineMultiDimension<SplineType> where SplineType : Spline {
        private readonly SplineType[] spline;

        /// <summary>コンストラクタ</summary>
        public SplineMultiDimension(int dimention, EndType type = EndType.Open) {
            if(dimention <= 0) {
                throw new ArgumentException();
            }

            this.spline = new SplineType[dimention];
            for(int i = 0; i < spline.Length; i++) {
                this.spline[i] = Spline.Construct<SplineType>(type);
            }
            this.Dimention = dimention;
        }

        /// <summary>制御点数</summary>
        public int Points => spline[0].Points;

        /// <summary>次元数</summary>
        public int Dimention { get; private set; }

        /// <summary>制御点を設定</summary>
        public void Set(params Vector[] v) {
            if(v == null) {
                throw new ArgumentNullException();
            }

            Set(v, v.Length);
        }

        /// <summary>制御点を設定</summary>
        public void Set(Vector[] v, int points) {
            Initialize();

            if(v == null) {
                throw new ArgumentNullException();
            }

            for(int i = 0; i < points; i++) {
                if(v[i].Dim != Dimention) {
                    throw new ArgumentException();
                }
            }

            double[] s = new double[points];

            for(int i = 0, j; i < Dimention; i++) {
                for(j = 0; j < points; j++) {
                    s[j] = v[j][i];
                }

                this.spline[i].Set(s);
            }
        }

        /// <summary>制御点を挿入</summary>
        public void Insert(int index, Vector new_v) {
            if(index < 0 || index > Points || new_v.Dim != Dimention) {
                throw new ArgumentException();
            }

            for(int i = 0; i < Dimention; i++) {
                this.spline[i].Insert(index, new_v[i]);
            }
        }

        /// <summary>制御点を除去</summary>
        public void Remove(int index) {
            if(index < 0 || index >= Points) {
                throw new ArgumentException();
            }

            for(int i = 0; i < Dimention; i++) {
                this.spline[i].Remove(index);
            }
        }

        /// <summary>初期化</summary>
        public void Initialize() {
            for(int i = 0; i < spline.Length; i++) {
                spline[i].Initialize();
            }
        }

        /// <summary>補間値</summary>
        public Vector Value(double t) {
            double[] v = new double[Dimention];

            for(int i = 0; i < Dimention; i++) {
                v[i] = spline[i].Value(t);
            }

            return new Vector(v);
        }

        /// <summary>補間微分値</summary>
        public Vector Diff(double t) {
            double[] v = new double[Dimention];

            for(int i = 0; i < Dimention; i++) {
                v[i] = spline[i].Diff(t);
            }

            return new Vector(v);
        }

        /// <summary>補間N次微分値</summary>
        public Vector Diff(double t, uint n) {
            double[] v = new double[Dimention];

            for(int i = 0; i < Dimention; i++) {
                v[i] = spline[i].Diff(t, n);
            }

            return new Vector(v);
        }

        /// <summary>等しいか判定</summary>
        public static bool operator ==(SplineMultiDimension<SplineType> s1, SplineMultiDimension<SplineType> s2) {
            if(ReferenceEquals(s1, s2)) {
                return true;
            }
            if((object)s1 == null || (object)s2 == null) {
                return false;
            }
            if(s1.Dimention != s2.Dimention) {
                return false;
            }
            for(int i = 0; i < s1.Dimention; i++) {
                if(s1.spline[i] != s2.spline[i]) {
                    return false;
                }
            }
            return true;
        }

        /// <summary>異なるか判定</summary>
        public static bool operator !=(SplineMultiDimension<SplineType> s1, SplineMultiDimension<SplineType> s2) {
            return !(s1 == s2);
        }

        /// <summary>等しいか判定</summary>
        public override bool Equals(object obj) {
            return obj is SplineMultiDimension<SplineType> ? (SplineMultiDimension<SplineType>)obj == this : false;
        }

        /// <summary>ハッシュ値</summary>
        public override int GetHashCode() {
            return spline[0].GetHashCode();
        }
    }
}