c# 二维 cubic 三次差值算法

二维 cubic 三次差值算法，非常好用，性能棒，

public  class Interpolate
    {
         /// <summary>
        /// 双三次差值算法。并行处理版。
        /// </summary>
        /// <param name="src">原数据</param>
        /// <param name="sizeSrc">元数据列数和行数</param>
        /// <param name="sizeDst">变换后列数和行数</param>
        /// <returns></returns>
        public  static double[] CubicInterpo(double[] src, int[] sizeSrc, int[] sizeDst)
        {
            int srcWidth = sizeSrc[0];
            int srcHeight = sizeSrc[1];

            int dstWidth = sizeDst[0];
            int dstHeight = sizeDst[1];

            // 计算比例系数
            double xFactor = (double)srcWidth / dstWidth;
            double yFactor = (double)srcHeight / dstHeight;

            double[] dstGrayData = new double[dstWidth * dstHeight];

            int ymax = srcHeight - 1;
            int xmax = srcWidth - 1;
            //int dst = 0;

            Parallel.For(0, dstHeight, y =>
            {
                // Y坐标
                double oy = (double)y * yFactor - 0.5;
                int oy1 = (int)oy;
                double dy = oy - (double)oy1;

                Parallel.For(0, dstWidth, x =>
                {
                    // X坐标
                    double ox = (double)x * xFactor - 0.5f;
                    int ox1 = (int)ox;
                    double dx = ox - (double)ox1;

                    // 像素值归零
                    double grayValue = 0;
                    for (int n = -1; n < 3; n++)
                    {
                        // Y系数
                        double k1 = BiCubicInterpolator(dy - (double)n);

                        int oy2 = oy1 + n;
                        if (oy2 < 0)
                        {
                            oy2 = 0;
                        }
                        if (oy2 > ymax)
                        {
                            oy2 = ymax;
                        }
                        for (int m = -1; m < 3; m++)
                        {
                            // X系数
                            double k2 = k1 * BiCubicInterpolator((double)m - dx);

                            int ox2 = ox1 + m;
                            if (ox2 < 0)
                            {
                                ox2 = 0;
                            }

                            if (ox2 > xmax)
                            {
                                ox2 = xmax;
                            }
                            grayValue += k2 * src[oy2 * srcWidth + ox2];
                        }
                    }

                    int dst = dstWidth * y + x;
                    dstGrayData[dst] = grayValue;
                });
            });
            return dstGrayData;
        }
        /// <summary>
        /// 双三次插值器。
        /// </summary>
        /// <param name="x">X Value.</param>
        /// <returns>Bicubic cooefficient.</returns>
        private static double BiCubicInterpolator(double x)
        {
            if (x < 0)
            {
                x = -x;
            }

            double biCoef = 0;

            if (x <= 1)
            {
                biCoef = (1.5 * x - 2.5) * x * x + 1;
            }
            else if (x < 2)
            {
                biCoef = ((-0.5 * x + 2.5) * x - 4) * x + 2;
            }

            return biCoef;
        }
    }