private byte[] BicubicResize(byte[] Data, int w, int h, double scalex, double scaley)
        {
            int w_new = (int)(w * scalex);
            int h_new = (int)(h * scaley);
            byte[] NewData = new byte[w_new * h_new * 4];
            byte c;
            int x, y;

            BicubicInterpolation interpolation = new BicubicInterpolation(scalex, scaley, w_new, h_new);
            for (c = 0; c < 4; c++)
                for (y = 0; y < h_new; y++)
                    for (x = 0; x < w_new; x++)
                        NewData[((x + y * w_new) << 2) | c] = interpolation.Bicubic(x, y, c, Data);
            return NewData;
        }

        public class BicubicInterpolation
        {
            double dx;
            double dy;
            double p0;
            double p1;
            double p2;
            double p3;
            double[] Rdx;
            double[] Rdy;
            double s;
            byte k;
            byte i;
            long pix0;
            long pix1;
            double[][] p;
            double scalex;
            double scaley;
            int w;
            int h;

            public BicubicInterpolation(double scalex = 1, double scaley = 1, int w = 0, int h = 0)
            {
                p = new double[4][];
                p[0] = new double[4];
                p[1] = new double[4];
                p[2] = new double[4];
                p[3] = new double[4];
                pix0 = 0;
                pix1 = 0;

                this.scalex = scalex;
                this.scaley = scaley;
                this.w = w;
                this.h = h;
            }

            public byte Bicubic(int x, int y, byte c, byte[] Data)
            {
                pix0 = (long)(x / scalex) + (long)(y / scaley) * (long)(w / scalex);
                for (k = 0; k < 4; k++)
                    for (i = 0; i < 4; i++)
                    {
                        pix1 = ((long)(pix0 + i + k * (w / scalex)) << 2) | c;
                        if (pix1 < Data.LongLength) p[i][k] = Data[pix1];
                        else p[i][k] = 0;
                    }

                //Temp = BicubicInterpolate(p, x - (long)(x / scalex), y - (long)(y / scaley));
                return (byte)Interpolation(x, y, p);
            }

            private double Interpolation(double tx, double ty, double[][] P)
            {
                dx = F(tx);
                dy = F(ty);

                Rdx = new double[4];

                Rdy = new double[4];
                for (int n = 0; n < 4; n++)
                {
                    Rdx[n] = R(n - 1 - dx);
                    Rdy[n] = R(n - 1 - dy);
                }

                s = 0;

                for (byte k = 0; k < 4; k++)
                    for (byte i = 0; i < 4; i++)
                        s += P[k][i] * Rdx[i] * Rdy[k];
                return s;
            }

            private double F(double x) => x - Math.Floor(x);
            private double P(double x) => (x > 0.0) ? x : 0.0;
            private double R(double x)
            {
                //return (Math.Pow(P(x + 2), 3) - 4 * Math.Pow(P(x + 1), 3) +
                //    6 * Math.Pow(P(x), 3) - 4 * Math.Pow(P(x - 1), 3)) / 6;
                p0 = P(x + 2);
                p1 = P(x + 1);
                p2 = P(x    );
                p3 = P(x - 1);
                return ((p0 * p0 * p0) - 4 * (p1 * p1 * p1) + 6 * (p2 * p2 * p2) - 4 * (p3 * p3 * p3)) / 6.0;
            }
        }