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];
            double[][] p = new double[4][];
            p[0] = new double[4];
            p[1] = new double[4];
            p[2] = new double[4];
            p[3] = new double[4];

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

        private double BicubicInterpolate(double[][] p, double x, double y)
        {
            double[] arr = new double[4];
            arr[0] = CubicInterpolate(p[0], y);
            arr[1] = CubicInterpolate(p[1], y);
            arr[2] = CubicInterpolate(p[2], y);
            arr[3] = CubicInterpolate(p[3], y);
            return CubicInterpolate(arr, x);
        }

        private double CubicInterpolate(double[] p, double x) => p[1] + 0.5 * x * (p[2] - p[0] + x *
            (2.0 * p[0] - 5.0 * p[1] + 4.0 * p[2] - p[3] + x * (3.0 * (p[1] - p[2]) + p[3] - p[0])));