跳到主要內容
黯羽輕揚每天積累一點點

js 實現曲線擬合(最小二乘法)

免費2015-06-12#JS#Solution#js曲线拟合#JavaScript曲线拟合#canvas平滑曲线

本文只給出可用的 js 代碼,對最小二乘法不作介紹

寫在前面

很多時候我們需要把一組點用平滑曲線連接,從點到曲線,就是曲線擬合。而最小二乘法是實現曲線擬合最簡單的方法,雖然效果不是很好,但勝在簡單

一。曲線擬合

用連續曲線近似地刻畫或擬擬平面上離散點組所表示的坐標之間的函數關係的一種數據處理方法

嗯,說白了就是用平滑曲線連接。手握鉛筆很容易搞定,但編程實現有一定難度,我們下面使用的方法叫最小二乘法,至於原理、數學公式、證明過程,抱歉,我不知道。

二。完整代碼

P.S. 下面的代碼是根據某 Java 實現改寫而來的,點我 查看原文

function PolyFitForcast() {
    /**
     * <p>
     * 函數功能:最小二乘法曲線擬合
     * </p>
     * <p>
     * 方程:Y = a(0) + a(1) * (X - X1)+ a(2) * (X - X1)^2 + ..... .+ a(m) * (X -
     * X1)^m X1 為 X 的平均數
     * </p>
     * 
     * @param x
     *            實型一維數組,長度為 n. 存放給定 n 個數據點的 X 坐標
     * @param y
     *            實型一維數組,長度為 n.存放給定 n 個數據點的 Y 坐標
     * @param n
     *            變量。給定數據點的個數
     * @param a
     *            實型一維數組,長度為 m.返回 m-1 次擬合多項式的 m 個係數
     * @param m
     *            擬合多項式的項數,即擬合多項式的最高次數為 m-1. 要求 m<=n 且 m<=20。若 m>n 或 m>20
     *            ,則本函數自動按 m=min{n,20} 處理.
     *            <p>
     *            Date:2007-12-25 16:21 PM
     *            </p>
     * @author qingbao-gao
     * @return 多項式係數存儲數組
     */
    function PolyFit(x, y, n, a, m) {
        var i, j, k;
        var z, p, c, g, q = 0, d1, d2;
        var s = new Array(20);
        var t = new Array(20);
        var b = new Array(20);
        var dt = new Array(3);
        for (i = 0; i <= m - 1; i++) {
            a[i] = 0.0;
        }
        if (m > n) {
            m = n;
        }
        if (m > 20) {
            m = 20;
        }
        z = 0.0;
        for (i = 0; i <= n - 1; i++) {
            z = z + x[i] / (1.0 * n);
        }
        b[0] = 1.0;
        d1 = 1.0 * n;
        p = 0.0;
        c = 0.0;
        for (i = 0; i <= n - 1; i++) {
            p = p + (x[i] - z);
            c = c + y[i];
        }
        c = c / d1;
        p = p / d1;
        a[0] = c * b[0];
        if (m > 1) {
            t[1] = 1.0;
            t[0] = -p;
            d2 = 0.0;
            c = 0.0;
            g = 0.0;
            for (i = 0; i <= n - 1; i++) {
                q = x[i] - z - p;
                d2 = d2 + q * q;
                c = c + y[i] * q;
                g = g + (x[i] - z) * q * q;
            }
            c = c / d2;
            p = g / d2;
            q = d2 / d1;
            d1 = d2;
            a[1] = c * t[1];
            a[0] = c * t[0] + a[0];
        }
        for (j = 2; j <= m - 1; j++) {
            s[j] = t[j - 1];
            s[j - 1] = -p * t[j - 1] + t[j - 2];
            if (j >= 3)
                for (k = j - 2; k >= 1; k--) {
                    s[k] = -p * t[k] + t[k - 1] - q * b[k];
                }
            s[0] = -p * t[0] - q * b[0];
            d2 = 0.0;
            c = 0.0;
            g = 0.0;
            for (i = 0; i <= n - 1; i++) {
                q = s[j];
                for (k = j - 1; k >= 0; k--) {
                    q = q * (x[i] - z) + s[k];
                }
                d2 = d2 + q * q;
                c = c + y[i] * q;
                g = g + (x[i] - z) * q * q;
            }
            c = c / d2;
            p = g / d2;
            q = d2 / d1;
            d1 = d2;
            a[j] = c * s[j];
            t[j] = s[j];
            for (k = j - 1; k >= 0; k--) {
                a[k] = c * s[k] + a[k];
                b[k] = t[k];
                t[k] = s[k];
            }
        }
        dt[0] = 0.0;
        dt[1] = 0.0;
        dt[2] = 0.0;
        for (i = 0; i <= n - 1; i++) {
            q = a[m - 1];
            for (k = m - 2; k >= 0; k--) {
                q = a[k] + q * (x[i] - z);
            }
            p = q - y[i];
            if (Math.abs(p) > dt[2]) {
                dt[2] = Math.abs(p);
            }
            dt[0] = dt[0] + p * p;
            dt[1] = dt[1] + Math.abs(p);
        }
        return a;
    }// end

    /**
     * <p>
     * 對 X 軸數據節點球平均值
     * </p>
     * 
     * @param x
     *            存儲 X 軸節點的數組
     *            <p>
     *            Date:2007-12-25 20:21 PM
     *            </p>
     * @author qingbao-gao
     * @return 平均值
     */
    function average(x) {
        var ave = 0;
        var sum = 0;
        if (x !== null) {
            for (var i = 0; i < x.length; i++) {
                sum += x[i];
            }
            ave = sum / x.length;
        }
        return ave;
    }

    /**
     * <p>
     * 由 X 值獲得 Y 值
     * </p>
     * @param x
     *            當前 X 軸輸入值,即為預測的月份
     * @param xx
     *            當前 X 軸輸入值的前 X 數據點
     * @param a
     *            存儲多項式係數的數組
     * @param m
     *            存儲多項式的最高次數的數組
     *            <p>
     *            Date:2007-12-25 PM 20:07
     *            </p>
     * @return 對應 X 軸節點值的 Y 軸值
     */
    function getY(x, xx, a, m) {
        var y = 0;
        var ave = average(xx);

        var l = 0;
        for (var i = 0; i < m; i++) {
            l = a[0];
            if (i > 0) {
                y += a[i] * Math.pow((x - ave), i);
            }
        }
        return (y + l);
    }

    /**
     * 返回擬合後的點坐標數組
     * @param  {Array} arr 點坐標數組
     * @return {Array}     擬合後的點坐標數組
     */
    this.get = function(arr) {
        var arrX = [], arrY = [];
        
        for (var i = 0; i < arr.length; i++) {
            arrX.push(arr[i].x);
            arrY.push(arr[i].y);
        }
        
        var len = arrY.length;
        var m = 3;
        var a = new Array(arrX.length);
        var aa = PolyFit(arrX, arrY, len, a, m);
        var arrRes = [];
        for(var i = 0; i < len; i++){
           arrRes.push({x: arrX[i], y: getY(arrX[i], arrX, aa, m)});
        }

        return arrRes;
    };
}

// var arr = [{x: 10, y: 10},{x: 40, y: 90},{x: 70, y: 65},{x: 100, y: 15}];
// // 最小二乘法擬合
// var res = new PolyFitForcast().get(arr);
// var canvas = document.createElement('canvas');
// var ctx2d = canvas.getContext('2d');

// ctx2d.lineWidth = 1;

// ctx2d.strokeStyle = '#000';
// ctx2d.beginPath();
// ctx2d.moveTo(arr[0].x, arr[0].y);
// ctx2d.lineTo(arr[1].x, arr[1].y);
// ctx2d.lineTo(arr[2].x, arr[2].y);
// ctx2d.lineTo(arr[3].x, arr[3].y);
// ctx2d.stroke();

// ctx2d.strokeStyle = '#f00';
// ctx2d.beginPath();
// ctx2d.moveTo(res[0].x, res[0].y);
// ctx2d.lineTo(res[1].x, res[1].y);
// ctx2d.lineTo(res[2].x, res[2].y);
// ctx2d.lineTo(res[3].x, res[3].y);
// ctx2d.stroke();
// document.body.appendChild(canvas);

三。效果

把上面代碼複製製到瀏覽器 console 去掉註釋再回車就能看到效果,當然,4 個點的效果很差,下面給出 10 個點擬合的效應圖:

curve

其實 10 個點擬合之後曲線還起來還不夠圓滑,但如果有一個場景恰好能夠提供很多個點,而且要求繪製速度還對精度要求不高的話,這種方式還是很不錯的

四。更平滑自然的曲線擬合

最小二乘法是沒指望了,除非反覆進行多次擬合,但效果提升並不明顯,曲線會越來越平,最後近似於直線,這不是我們想要的結果

還有一種方法叫貝塞爾曲線,除了一組信息點,還需要提供兩個控制點,如果控制點選取得當就能實現相當完美的平滑曲線效果,難點是如何根據現有的信息點求出 2 個控制點,具體方法請查看 Smoother Signatures

評論

暫無評論,快來發表你的看法吧

提交評論