C# spline2dinterpolant类代码示例

        /*************************************************************************
        This subroutine performs linear transformation of the spline.

        Input parameters:
            C   -   spline interpolant.
            A, B-   transformation coefficients: S2(x,y) = A*S(x,y) + B
            
        Output parameters:
            C   -   transformed spline

          -- ALGLIB PROJECT --
             Copyright 30.06.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dlintransf(spline2dinterpolant c,
            double a,
            double b)
        {
            double[] x = new double[0];
            double[] y = new double[0];
            double[] f = new double[0];
            int i = 0;
            int j = 0;

            alglib.ap.assert(c.stype==-3 || c.stype==-1, "Spline2DLinTransF: incorrect C (incorrect parameter C.SType)");
            x = new double[c.n];
            y = new double[c.m];
            f = new double[c.m*c.n*c.d];
            for(j=0; j<=c.n-1; j++)
            {
                x[j] = c.x[j];
            }
            for(i=0; i<=c.m-1; i++)
            {
                y[i] = c.y[i];
            }
            for(i=0; i<=c.m*c.n*c.d-1; i++)
            {
                f[i] = a*c.f[i]+b;
            }
            if( c.stype==-3 )
            {
                spline2dbuildbicubicv(x, c.n, y, c.m, f, c.d, c);
            }
            if( c.stype==-1 )
            {
                spline2dbuildbilinearv(x, c.n, y, c.m, f, c.d, c);
            }
        }

        /*************************************************************************
        This subroutine calculates the value of the bilinear or bicubic spline  at
        the given point X and its derivatives.

        Input parameters:
            C   -   spline interpolant.
            X, Y-   point

        Output parameters:
            F   -   S(x,y)
            FX  -   dS(x,y)/dX
            FY  -   dS(x,y)/dY
            FXY -   d2S(x,y)/dXdY

          -- ALGLIB PROJECT --
             Copyright 05.07.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2ddiff(spline2dinterpolant c,
            double x,
            double y,
            ref double f,
            ref double fx,
            ref double fy,
            ref double fxy)
        {
            double t = 0;
            double dt = 0;
            double u = 0;
            double du = 0;
            int ix = 0;
            int iy = 0;
            int l = 0;
            int r = 0;
            int h = 0;
            int s1 = 0;
            int s2 = 0;
            int s3 = 0;
            int s4 = 0;
            int sfx = 0;
            int sfy = 0;
            int sfxy = 0;
            double y1 = 0;
            double y2 = 0;
            double y3 = 0;
            double y4 = 0;
            double v = 0;
            double t0 = 0;
            double t1 = 0;
            double t2 = 0;
            double t3 = 0;
            double u0 = 0;
            double u1 = 0;
            double u2 = 0;
            double u3 = 0;

            f = 0;
            fx = 0;
            fy = 0;
            fxy = 0;

            alglib.ap.assert(c.stype==-1 || c.stype==-3, "Spline2DDiff: incorrect C (incorrect parameter C.SType)");
            alglib.ap.assert(math.isfinite(x) && math.isfinite(y), "Spline2DDiff: X or Y contains NaN or Infinite value");
            
            //
            // Prepare F, dF/dX, dF/dY, d2F/dXdY
            //
            f = 0;
            fx = 0;
            fy = 0;
            fxy = 0;
            if( c.d!=1 )
            {
                return;
            }
            
            //
            // Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
            //
            l = 0;
            r = c.n-1;
            while( l!=r-1 )
            {
                h = (l+r)/2;
                if( (double)(c.x[h])>=(double)(x) )
                {
                    r = h;
                }
                else
                {
                    l = h;
                }
            }
            t = (x-c.x[l])/(c.x[l+1]-c.x[l]);
            dt = 1.0/(c.x[l+1]-c.x[l]);
            ix = l;
            
            //
            // Binary search in the [ y[0], ..., y[m-2] ] (y[m-1] is not included)
            //
            l = 0;
//.........这里部分代码省略.........

        /*************************************************************************
        This subroutine performs linear transformation of the spline argument.

        Input parameters:
            C       -   spline interpolant
            AX, BX  -   transformation coefficients: x = A*t + B
            AY, BY  -   transformation coefficients: y = A*u + B
        Result:
            C   -   transformed spline

          -- ALGLIB PROJECT --
             Copyright 30.06.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dlintransxy(spline2dinterpolant c,
            double ax,
            double bx,
            double ay,
            double by)
        {
            double[] x = new double[0];
            double[] y = new double[0];
            double[] f = new double[0];
            double[] v = new double[0];
            int i = 0;
            int j = 0;
            int k = 0;

            alglib.ap.assert(c.stype==-3 || c.stype==-1, "Spline2DLinTransXY: incorrect C (incorrect parameter C.SType)");
            alglib.ap.assert(math.isfinite(ax), "Spline2DLinTransXY: AX is infinite or NaN");
            alglib.ap.assert(math.isfinite(bx), "Spline2DLinTransXY: BX is infinite or NaN");
            alglib.ap.assert(math.isfinite(ay), "Spline2DLinTransXY: AY is infinite or NaN");
            alglib.ap.assert(math.isfinite(by), "Spline2DLinTransXY: BY is infinite or NaN");
            x = new double[c.n];
            y = new double[c.m];
            f = new double[c.m*c.n*c.d];
            for(j=0; j<=c.n-1; j++)
            {
                x[j] = c.x[j];
            }
            for(i=0; i<=c.m-1; i++)
            {
                y[i] = c.y[i];
            }
            for(i=0; i<=c.m-1; i++)
            {
                for(j=0; j<=c.n-1; j++)
                {
                    for(k=0; k<=c.d-1; k++)
                    {
                        f[c.d*(i*c.n+j)+k] = c.f[c.d*(i*c.n+j)+k];
                    }
                }
            }
            
            //
            // Handle different combinations of AX/AY
            //
            if( (double)(ax)==(double)(0) && (double)(ay)!=(double)(0) )
            {
                for(i=0; i<=c.m-1; i++)
                {
                    spline2dcalcvbuf(c, bx, y[i], ref v);
                    y[i] = (y[i]-by)/ay;
                    for(j=0; j<=c.n-1; j++)
                    {
                        for(k=0; k<=c.d-1; k++)
                        {
                            f[c.d*(i*c.n+j)+k] = v[k];
                        }
                    }
                }
            }
            if( (double)(ax)!=(double)(0) && (double)(ay)==(double)(0) )
            {
                for(j=0; j<=c.n-1; j++)
                {
                    spline2dcalcvbuf(c, x[j], by, ref v);
                    x[j] = (x[j]-bx)/ax;
                    for(i=0; i<=c.m-1; i++)
                    {
                        for(k=0; k<=c.d-1; k++)
                        {
                            f[c.d*(i*c.n+j)+k] = v[k];
                        }
                    }
                }
            }
            if( (double)(ax)!=(double)(0) && (double)(ay)!=(double)(0) )
            {
                for(j=0; j<=c.n-1; j++)
                {
                    x[j] = (x[j]-bx)/ax;
                }
                for(i=0; i<=c.m-1; i++)
                {
                    y[i] = (y[i]-by)/ay;
                }
            }
            if( (double)(ax)==(double)(0) && (double)(ay)==(double)(0) )
            {
//.........这里部分代码省略.........

        /*************************************************************************
        This subroutine was deprecated in ALGLIB 3.6.0

        We recommend you to switch  to  Spline2DBuildBicubicV(),  which  is  more
        flexible and accepts its arguments in more convenient order.

          -- ALGLIB PROJECT --
             Copyright 05.07.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dbuildbicubic(double[] x,
            double[] y,
            double[,] f,
            int m,
            int n,
            spline2dinterpolant c)
        {
            int sfx = 0;
            int sfy = 0;
            int sfxy = 0;
            double[,] dx = new double[0,0];
            double[,] dy = new double[0,0];
            double[,] dxy = new double[0,0];
            double t = 0;
            int i = 0;
            int j = 0;
            int k = 0;

            f = (double[,])f.Clone();

            alglib.ap.assert(n>=2, "Spline2DBuildBicubicSpline: N<2");
            alglib.ap.assert(m>=2, "Spline2DBuildBicubicSpline: M<2");
            alglib.ap.assert(alglib.ap.len(x)>=n && alglib.ap.len(y)>=m, "Spline2DBuildBicubic: length of X or Y is too short (Length(X/Y)<N/M)");
            alglib.ap.assert(apserv.isfinitevector(x, n) && apserv.isfinitevector(y, m), "Spline2DBuildBicubic: X or Y contains NaN or Infinite value");
            alglib.ap.assert(alglib.ap.rows(f)>=m && alglib.ap.cols(f)>=n, "Spline2DBuildBicubic: size of F is too small (rows(F)<M or cols(F)<N)");
            alglib.ap.assert(apserv.apservisfinitematrix(f, m, n), "Spline2DBuildBicubic: F contains NaN or Infinite value");
            
            //
            // Fill interpolant:
            //  F[0]...F[N*M-1]:
            //      f(i,j) table. f(0,0), f(0, 1), f(0,2) and so on...
            //  F[N*M]...F[2*N*M-1]:
            //      df(i,j)/dx table.
            //  F[2*N*M]...F[3*N*M-1]:
            //      df(i,j)/dy table.
            //  F[3*N*M]...F[4*N*M-1]:
            //      d2f(i,j)/dxdy table.
            //
            c.k = 3;
            c.d = 1;
            c.n = n;
            c.m = m;
            c.stype = -3;
            sfx = c.n*c.m;
            sfy = 2*c.n*c.m;
            sfxy = 3*c.n*c.m;
            c.x = new double[c.n];
            c.y = new double[c.m];
            c.f = new double[4*c.n*c.m];
            for(i=0; i<=c.n-1; i++)
            {
                c.x[i] = x[i];
            }
            for(i=0; i<=c.m-1; i++)
            {
                c.y[i] = y[i];
            }
            
            //
            // Sort points
            //
            for(j=0; j<=c.n-1; j++)
            {
                k = j;
                for(i=j+1; i<=c.n-1; i++)
                {
                    if( (double)(c.x[i])<(double)(c.x[k]) )
                    {
                        k = i;
                    }
                }
                if( k!=j )
                {
                    for(i=0; i<=c.m-1; i++)
                    {
                        t = f[i,j];
                        f[i,j] = f[i,k];
                        f[i,k] = t;
                    }
                    t = c.x[j];
                    c.x[j] = c.x[k];
                    c.x[k] = t;
                }
            }
            for(i=0; i<=c.m-1; i++)
            {
                k = i;
                for(j=i+1; j<=c.m-1; j++)
                {
                    if( (double)(c.y[j])<(double)(c.y[k]) )
                    {
//.........这里部分代码省略.........

        /*************************************************************************
        This subroutine makes the copy of the spline model.

        Input parameters:
            C   -   spline interpolant

        Output parameters:
            CC  -   spline copy

          -- ALGLIB PROJECT --
             Copyright 29.06.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dcopy(ref spline2dinterpolant c,
            ref spline2dinterpolant cc)
        {
            int n = 0;
            int i_ = 0;

            System.Diagnostics.Debug.Assert(c.k==1 | c.k==3, "Spline2DCopy: incorrect C!");
            cc.k = c.k;
            n = (int)Math.Round(c.c[0]);
            cc.c = new double[n];
            for(i_=0; i_<=n-1;i_++)
            {
                cc.c[i_] = c.c[i_];
            }
        }

        /*************************************************************************
        This subroutine builds bicubic vector-valued spline.

        Input parameters:
            X   -   spline abscissas, array[0..N-1]
            Y   -   spline ordinates, array[0..M-1]
            F   -   function values, array[0..M*N*D-1]:
                    * first D elements store D values at (X[0],Y[0])
                    * next D elements store D values at (X[1],Y[0])
                    * general form - D function values at (X[i],Y[j]) are stored
                      at F[D*(J*N+I)...D*(J*N+I)+D-1].
            M,N -   grid size, M>=2, N>=2
            D   -   vector dimension, D>=1

        Output parameters:
            C   -   spline interpolant

          -- ALGLIB PROJECT --
             Copyright 16.04.2012 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dbuildbicubicv(double[] x,
            int n,
            double[] y,
            int m,
            double[] f,
            int d,
            spline2dinterpolant c)
        {
            double[,] tf = new double[0,0];
            double[,] dx = new double[0,0];
            double[,] dy = new double[0,0];
            double[,] dxy = new double[0,0];
            double t = 0;
            int i = 0;
            int j = 0;
            int k = 0;
            int di = 0;

            f = (double[])f.Clone();

            alglib.ap.assert(n>=2, "Spline2DBuildBicubicV: N is less than 2");
            alglib.ap.assert(m>=2, "Spline2DBuildBicubicV: M is less than 2");
            alglib.ap.assert(d>=1, "Spline2DBuildBicubicV: invalid argument D (D<1)");
            alglib.ap.assert(alglib.ap.len(x)>=n && alglib.ap.len(y)>=m, "Spline2DBuildBicubicV: length of X or Y is too short (Length(X/Y)<N/M)");
            alglib.ap.assert(apserv.isfinitevector(x, n) && apserv.isfinitevector(y, m), "Spline2DBuildBicubicV: X or Y contains NaN or Infinite value");
            k = n*m*d;
            alglib.ap.assert(alglib.ap.len(f)>=k, "Spline2DBuildBicubicV: length of F is too short (Length(F)<N*M*D)");
            alglib.ap.assert(apserv.isfinitevector(f, k), "Spline2DBuildBicubicV: F contains NaN or Infinite value");
            
            //
            // Fill interpolant:
            //  F[0]...F[N*M*D-1]:
            //      f(i,j) table. f(0,0), f(0, 1), f(0,2) and so on...
            //  F[N*M*D]...F[2*N*M*D-1]:
            //      df(i,j)/dx table.
            //  F[2*N*M*D]...F[3*N*M*D-1]:
            //      df(i,j)/dy table.
            //  F[3*N*M*D]...F[4*N*M*D-1]:
            //      d2f(i,j)/dxdy table.
            //
            c.k = 3;
            c.d = d;
            c.n = n;
            c.m = m;
            c.stype = -3;
            k = 4*k;
            c.x = new double[c.n];
            c.y = new double[c.m];
            c.f = new double[k];
            tf = new double[c.m, c.n];
            for(i=0; i<=c.n-1; i++)
            {
                c.x[i] = x[i];
            }
            for(i=0; i<=c.m-1; i++)
            {
                c.y[i] = y[i];
            }
            
            //
            // Sort points
            //
            for(j=0; j<=c.n-1; j++)
            {
                k = j;
                for(i=j+1; i<=c.n-1; i++)
                {
                    if( (double)(c.x[i])<(double)(c.x[k]) )
                    {
                        k = i;
                    }
                }
                if( k!=j )
                {
                    for(i=0; i<=c.m-1; i++)
                    {
                        for(di=0; di<=c.d-1; di++)
                        {
                            t = f[c.d*(i*c.n+j)+di];
                            f[c.d*(i*c.n+j)+di] = f[c.d*(i*c.n+k)+di];
//.........这里部分代码省略.........

        /*************************************************************************
        This subroutine calculates bilinear or bicubic vector-valued spline at the
        given point (X,Y).

        INPUT PARAMETERS:
            C   -   spline interpolant.
            X, Y-   point

        OUTPUT PARAMETERS:
            F   -   array[D] which stores function values.  F is out-parameter and
                    it  is  reallocated  after  call to this function. In case you
                    want  to    reuse  previously  allocated  F,   you   may   use
                    Spline2DCalcVBuf(),  which  reallocates  F only when it is too
                    small.

          -- ALGLIB PROJECT --
             Copyright 16.04.2012 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dcalcv(spline2dinterpolant c,
            double x,
            double y,
            ref double[] f)
        {
            f = new double[0];

            alglib.ap.assert(c.stype==-1 || c.stype==-3, "Spline2DCalcV: incorrect C (incorrect parameter C.SType)");
            alglib.ap.assert(math.isfinite(x) && math.isfinite(y), "Spline2DCalcV: either X=NaN/Infinite or Y=NaN/Infinite");
            f = new double[c.d];
            spline2dcalcvbuf(c, x, y, ref f);
        }

            /*************************************************************************
            This subroutine unpacks two-dimensional spline into the coefficients table

            Input parameters:
                C   -   spline interpolant.

            Result:
                M, N-   grid size (x-axis and y-axis)
                Tbl -   coefficients table, unpacked format,
                        [0..(N-1)*(M-1)-1, 0..19].
                        For I = 0...M-2, J=0..N-2:
                            K =  I*(N-1)+J
                            Tbl[K,0] = X[j]
                            Tbl[K,1] = X[j+1]
                            Tbl[K,2] = Y[i]
                            Tbl[K,3] = Y[i+1]
                            Tbl[K,4] = C00
                            Tbl[K,5] = C01
                            Tbl[K,6] = C02
                            Tbl[K,7] = C03
                            Tbl[K,8] = C10
                            Tbl[K,9] = C11
                            ...
                            Tbl[K,19] = C33
                        On each grid square spline is equals to:
                            S(x) = SUM(c[i,j]*(x^i)*(y^j), i=0..3, j=0..3)
                            t = x-x[j]
                            u = y-y[i]

              -- ALGLIB PROJECT --
                 Copyright 29.06.2007 by Bochkanov Sergey
            *************************************************************************/
            public static void spline2dunpack(spline2dinterpolant c,
                ref int m,
                ref int n,
                ref double[,] tbl) {
                int i = 0;
                int j = 0;
                int ci = 0;
                int cj = 0;
                int k = 0;
                int p = 0;
                int shift = 0;
                int s1 = 0;
                int s2 = 0;
                int s3 = 0;
                int s4 = 0;
                int sf = 0;
                int sfx = 0;
                int sfy = 0;
                int sfxy = 0;
                double y1 = 0;
                double y2 = 0;
                double y3 = 0;
                double y4 = 0;
                double dt = 0;
                double du = 0;

                m = 0;
                n = 0;
                tbl = new double[0, 0];

                ap.assert((int)Math.Round(c.c[1]) == -3 | (int)Math.Round(c.c[1]) == -1, "SplineUnpack2D: incorrect C!");
                n = (int)Math.Round(c.c[2]);
                m = (int)Math.Round(c.c[3]);
                tbl = new double[(n - 1) * (m - 1) - 1 + 1, 19 + 1];

                //
                // Fill
                //
                for(i = 0; i <= m - 2; i++) {
                    for(j = 0; j <= n - 2; j++) {
                        p = i * (n - 1) + j;
                        tbl[p, 0] = c.c[4 + j];
                        tbl[p, 1] = c.c[4 + j + 1];
                        tbl[p, 2] = c.c[4 + n + i];
                        tbl[p, 3] = c.c[4 + n + i + 1];
                        dt = 1 / (tbl[p, 1] - tbl[p, 0]);
                        du = 1 / (tbl[p, 3] - tbl[p, 2]);

                        //
                        // Bilinear interpolation
                        //
                        if((int)Math.Round(c.c[1]) == -1) {
                            for(k = 4; k <= 19; k++) {
                                tbl[p, k] = 0;
                            }
                            shift = 4 + n + m;
                            y1 = c.c[shift + n * i + j];
                            y2 = c.c[shift + n * i + (j + 1)];
                            y3 = c.c[shift + n * (i + 1) + (j + 1)];
                            y4 = c.c[shift + n * (i + 1) + j];
                            tbl[p, 4] = y1;
                            tbl[p, 4 + 1 * 4 + 0] = y2 - y1;
                            tbl[p, 4 + 0 * 4 + 1] = y4 - y1;
                            tbl[p, 4 + 1 * 4 + 1] = y3 - y2 - y4 + y1;
                        }

                        //
                        // Bicubic interpolation
//.........这里部分代码省略.........

            /*************************************************************************
            This subroutine performs linear transformation of the spline argument.

            Input parameters:
                C       -   spline interpolant
                AX, BX  -   transformation coefficients: x = A*t + B
                AY, BY  -   transformation coefficients: y = A*u + B
            Result:
                C   -   transformed spline

              -- ALGLIB PROJECT --
                 Copyright 30.06.2007 by Bochkanov Sergey
            *************************************************************************/
            public static void spline2dlintransxy(spline2dinterpolant c,
                double ax,
                double bx,
                double ay,
                double by) {
                int i = 0;
                int j = 0;
                int n = 0;
                int m = 0;
                double v = 0;
                double[] x = new double[0];
                double[] y = new double[0];
                double[,] f = new double[0, 0];
                int typec = 0;

                typec = (int)Math.Round(c.c[1]);
                ap.assert(typec == -3 | typec == -1, "Spline2DLinTransXY: incorrect C!");
                n = (int)Math.Round(c.c[2]);
                m = (int)Math.Round(c.c[3]);
                x = new double[n - 1 + 1];
                y = new double[m - 1 + 1];
                f = new double[m - 1 + 1, n - 1 + 1];
                for(j = 0; j <= n - 1; j++) {
                    x[j] = c.c[4 + j];
                }
                for(i = 0; i <= m - 1; i++) {
                    y[i] = c.c[4 + n + i];
                }
                for(i = 0; i <= m - 1; i++) {
                    for(j = 0; j <= n - 1; j++) {
                        f[i, j] = c.c[4 + n + m + i * n + j];
                    }
                }

                //
                // Special case: AX=0 or AY=0
                //
                if((double)(ax) == (double)(0)) {
                    for(i = 0; i <= m - 1; i++) {
                        v = spline2dcalc(c, bx, y[i]);
                        for(j = 0; j <= n - 1; j++) {
                            f[i, j] = v;
                        }
                    }
                    if(typec == -3) {
                        spline2dbuildbicubic(x, y, f, m, n, c);
                    }
                    if(typec == -1) {
                        spline2dbuildbilinear(x, y, f, m, n, c);
                    }
                    ax = 1;
                    bx = 0;
                }
                if((double)(ay) == (double)(0)) {
                    for(j = 0; j <= n - 1; j++) {
                        v = spline2dcalc(c, x[j], by);
                        for(i = 0; i <= m - 1; i++) {
                            f[i, j] = v;
                        }
                    }
                    if(typec == -3) {
                        spline2dbuildbicubic(x, y, f, m, n, c);
                    }
                    if(typec == -1) {
                        spline2dbuildbilinear(x, y, f, m, n, c);
                    }
                    ay = 1;
                    by = 0;
                }

                //
                // General case: AX<>0, AY<>0
                // Unpack, scale and pack again.
                //
                for(j = 0; j <= n - 1; j++) {
                    x[j] = (x[j] - bx) / ax;
                }
                for(i = 0; i <= m - 1; i++) {
                    y[i] = (y[i] - by) / ay;
                }
                if(typec == -3) {
                    spline2dbuildbicubic(x, y, f, m, n, c);
                }
                if(typec == -1) {
                    spline2dbuildbilinear(x, y, f, m, n, c);
                }
            }

            /*************************************************************************
            This subroutine calculates the value of the bilinear or bicubic spline  at
            the given point X.

            Input parameters:
                C   -   coefficients table.
                        Built by BuildBilinearSpline or BuildBicubicSpline.
                X, Y-   point

            Result:
                S(x,y)

              -- ALGLIB PROJECT --
                 Copyright 05.07.2007 by Bochkanov Sergey
            *************************************************************************/
            public static double spline2dcalc(spline2dinterpolant c,
                double x,
                double y) {
                double result = 0;
                double v = 0;
                double vx = 0;
                double vy = 0;
                double vxy = 0;

                spline2ddiff(c, x, y, ref v, ref vx, ref vy, ref vxy);
                result = v;
                return result;
            }

            /*************************************************************************
            This subroutine builds bicubic spline coefficients table.

            Input parameters:
                X   -   spline abscissas, array[0..N-1]
                Y   -   spline ordinates, array[0..M-1]
                F   -   function values, array[0..M-1,0..N-1]
                M,N -   grid size, M>=2, N>=2

            Output parameters:
                C   -   spline interpolant

              -- ALGLIB PROJECT --
                 Copyright 05.07.2007 by Bochkanov Sergey
            *************************************************************************/
            public static void spline2dbuildbicubic(double[] x,
                double[] y,
                double[,] f,
                int m,
                int n,
                spline2dinterpolant c) {
                int i = 0;
                int j = 0;
                int k = 0;
                int tblsize = 0;
                int shift = 0;
                double t = 0;
                double[,] dx = new double[0, 0];
                double[,] dy = new double[0, 0];
                double[,] dxy = new double[0, 0];

                x = (double[])x.Clone();
                y = (double[])y.Clone();
                f = (double[,])f.Clone();

                ap.assert(n >= 2 & m >= 2, "BuildBicubicSpline: N<2 or M<2!");

                //
                // Sort points
                //
                for(j = 0; j <= n - 1; j++) {
                    k = j;
                    for(i = j + 1; i <= n - 1; i++) {
                        if((double)(x[i]) < (double)(x[k])) {
                            k = i;
                        }
                    }
                    if(k != j) {
                        for(i = 0; i <= m - 1; i++) {
                            t = f[i, j];
                            f[i, j] = f[i, k];
                            f[i, k] = t;
                        }
                        t = x[j];
                        x[j] = x[k];
                        x[k] = t;
                    }
                }
                for(i = 0; i <= m - 1; i++) {
                    k = i;
                    for(j = i + 1; j <= m - 1; j++) {
                        if((double)(y[j]) < (double)(y[k])) {
                            k = j;
                        }
                    }
                    if(k != i) {
                        for(j = 0; j <= n - 1; j++) {
                            t = f[i, j];
                            f[i, j] = f[k, j];
                            f[k, j] = t;
                        }
                        t = y[i];
                        y[i] = y[k];
                        y[k] = t;
                    }
                }

                //
                // Fill C:
                //  C[0]            -   length(C)
                //  C[1]            -   type(C):
                //                      -1 = bilinear interpolant
                //                           (see BuildBilinearInterpolant)
                //                      -3 = general cubic spline
                //  C[2]:
                //      N (x count)
                //  C[3]:
                //      M (y count)
                //  C[4]...C[4+N-1]:
                //      x[i], i = 0...N-1
                //  C[4+N]...C[4+N+M-1]:
                //      y[i], i = 0...M-1
                //  C[4+N+M]...C[4+N+M+(N*M-1)]:
                //      f(i,j) table. f(0,0), f(0, 1), f(0,2) and so on...
                //  C[4+N+M+N*M]...C[4+N+M+(2*N*M-1)]:
                //      df(i,j)/dx table.
                //  C[4+N+M+2*N*M]...C[4+N+M+(3*N*M-1)]:
                //      df(i,j)/dy table.
                //  C[4+N+M+3*N*M]...C[4+N+M+(4*N*M-1)]:
                //      d2f(i,j)/dxdy table.
//.........这里部分代码省略.........

        /*************************************************************************
        Unserialization of the spline interpolant

        INPUT PARAMETERS:
            RA  -   array of real numbers which contains interpolant,

        OUTPUT PARAMETERS:
            B   -   spline interpolant

          -- ALGLIB --
             Copyright 17.08.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dunserialize(ref double[] ra,
            ref spline2dinterpolant c)
        {
            int clen = 0;
            int i_ = 0;
            int i1_ = 0;

            System.Diagnostics.Debug.Assert((int)Math.Round(ra[1])==spline2dvnum, "Spline2DUnserialize: corrupted array!");
            c.k = (int)Math.Round(ra[2]);
            clen = (int)Math.Round(ra[3]);
            c.c = new double[clen];
            i1_ = (3) - (0);
            for(i_=0; i_<=clen-1;i_++)
            {
                c.c[i_] = ra[i_+i1_];
            }
        }

        /*************************************************************************
        Serialization of the spline interpolant

        INPUT PARAMETERS:
            B   -   spline interpolant

        OUTPUT PARAMETERS:
            RA      -   array of real numbers which contains interpolant,
                        array[0..RLen-1]
            RLen    -   RA lenght

          -- ALGLIB --
             Copyright 17.08.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dserialize(ref spline2dinterpolant c,
            ref double[] ra,
            ref int ralen)
        {
            int clen = 0;
            int i_ = 0;
            int i1_ = 0;

            System.Diagnostics.Debug.Assert(c.k==1 | c.k==3, "Spline2DSerialize: incorrect C!");
            clen = (int)Math.Round(c.c[0]);
            ralen = 3+clen;
            ra = new double[ralen];
            ra[0] = ralen;
            ra[1] = spline2dvnum;
            ra[2] = c.k;
            i1_ = (0) - (3);
            for(i_=3; i_<=3+clen-1;i_++)
            {
                ra[i_] = c.c[i_+i1_];
            }
        }

        /*************************************************************************
        This subroutine makes the copy of the spline model.

        Input parameters:
            C   -   spline interpolant

        Output parameters:
            CC  -   spline copy

          -- ALGLIB PROJECT --
             Copyright 29.06.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dcopy(spline2dinterpolant c,
            spline2dinterpolant cc)
        {
            int tblsize = 0;
            int i_ = 0;

            alglib.ap.assert(c.k==1 || c.k==3, "Spline2DCopy: incorrect C (incorrect parameter C.K)");
            cc.k = c.k;
            cc.n = c.n;
            cc.m = c.m;
            cc.d = c.d;
            cc.stype = c.stype;
            tblsize = -1;
            if( c.stype==-3 )
            {
                tblsize = 4*c.n*c.m*c.d;
            }
            if( c.stype==-1 )
            {
                tblsize = c.n*c.m*c.d;
            }
            alglib.ap.assert(tblsize>0, "Spline2DCopy: internal error");
            cc.x = new double[cc.n];
            cc.y = new double[cc.m];
            cc.f = new double[tblsize];
            for(i_=0; i_<=cc.n-1;i_++)
            {
                cc.x[i_] = c.x[i_];
            }
            for(i_=0; i_<=cc.m-1;i_++)
            {
                cc.y[i_] = c.y[i_];
            }
            for(i_=0; i_<=tblsize-1;i_++)
            {
                cc.f[i_] = c.f[i_];
            }
        }

            /*************************************************************************
            This subroutine performs linear transformation of the spline.

            Input parameters:
                C   -   spline interpolant.
                A, B-   transformation coefficients: S2(x,y) = A*S(x,y) + B
            
            Output parameters:
                C   -   transformed spline

              -- ALGLIB PROJECT --
                 Copyright 30.06.2007 by Bochkanov Sergey
            *************************************************************************/
            public static void spline2dlintransf(spline2dinterpolant c,
                double a,
                double b) {
                int i = 0;
                int j = 0;
                int n = 0;
                int m = 0;
                double[] x = new double[0];
                double[] y = new double[0];
                double[,] f = new double[0, 0];
                int typec = 0;

                typec = (int)Math.Round(c.c[1]);
                ap.assert(typec == -3 | typec == -1, "Spline2DLinTransXY: incorrect C!");
                n = (int)Math.Round(c.c[2]);
                m = (int)Math.Round(c.c[3]);
                x = new double[n - 1 + 1];
                y = new double[m - 1 + 1];
                f = new double[m - 1 + 1, n - 1 + 1];
                for(j = 0; j <= n - 1; j++) {
                    x[j] = c.c[4 + j];
                }
                for(i = 0; i <= m - 1; i++) {
                    y[i] = c.c[4 + n + i];
                }
                for(i = 0; i <= m - 1; i++) {
                    for(j = 0; j <= n - 1; j++) {
                        f[i, j] = a * c.c[4 + n + m + i * n + j] + b;
                    }
                }
                if(typec == -3) {
                    spline2dbuildbicubic(x, y, f, m, n, c);
                }
                if(typec == -1) {
                    spline2dbuildbilinear(x, y, f, m, n, c);
                }
            }