Evaluating cardinal splines

cspline.js

/*global sampler, drawPoints, pointString, lastOf */

function CSpline(points, variance) {
    this.tension = 0.2;
    this.points = points;
    this.variance = variance || points.map(function () { return 0; });
    this.interval = this.points[this._nextpi(0)].x - this.points[0].x;
    this.tangent = points.map(function (d, i) {
        var previ = this._prevpi(i),
            nexti = this._nextpi(i);
        //var previ = i - 1 < 0 ? 0 : i - 1,
        //    nexti = i + 1 >= this.points.length ? this.points.length - 1 : i + 1;
        return [
            (1 - this.tension) *
                (this.points[nexti].x - this.points[previ].x),
            (1 - this.tension) *
                //(this._getvalue(nexti) - this._getvalue(previ))
                (this.points[nexti].y - this.points[previ].y)
        ];
    }, this);
}

CSpline.prototype._getvalue = function (i) {
    return this.points[i].y - this.variance[i];
};

CSpline.prototype._prevpi = function (i) {
    var p = i - 1;
    while (p >= 0 && Math.abs(this.points[i].x - this.points[p].x) < 1e-6) {
        p--;
    }
    if (p < 0) {
        return i;
    }
    return p;
};

CSpline.prototype._nextpi = function (i) {
    var p = i + 1,
        N = this.points.length;
    while (p < N && Math.abs(this.points[p].x - this.points[i].x) < 1e-6) {
        p++;
    }
    if (p >= N) {
        return i;
    }
    return p;
};

CSpline.prototype._h00 = function (t) {
    return 2 * t * t * t - 3 * t * t + 1;
};

CSpline.prototype._h10 = function (t) {
    return t * t * t - 2 * t * t + t;
};

CSpline.prototype._h01 = function (t) {
    return -2 * t * t * t + 3 * t * t;
};

CSpline.prototype._h11 = function (t) {
    return t * t * t - t * t;
};

CSpline.prototype._findIndex = function (x) {
    var left = 0,
        right = this.points.length,
        mid;
    while (left < right) {
        mid = ~~((left + right) / 2);
        if (this.points[mid].x <= x &&
            (mid === this.points.length - 1 ||
                this.points[mid + 1].x >= x)) {
            return mid;
        } else if (x < this.points[mid].x) {
            right = mid;
        } else {
            left = mid + 1;
        }
    }
    return left;
};

CSpline.prototype.evaluate = function (x) {
    var i = this._findIndex(x);
    if (i >= this.points.length - 1) {
        i = this._prevpi(this.points.length - 1);
    }
    var npi = this._nextpi(i),
        x0 = this.points[i].x,
        x1 = this.points[npi].x,
        d0 = this.tangent[i][0],
        d1 = this.tangent[npi][0],
        a = 2 * x0 + d0 - 2 * x1 + d1,
        b = -3 * x0 - 2 * d0 + 3 * x1 - d1,
        c = d0,
        d = x0 - x,
        val = function (m) {
            return a * m * m * m + b * m * m + c * m + d;
        },
        sig = function (m) {
            return m > 0 ? 1 : -1;
        },
        l = 0, lv = val(0),
        r = 1, rv = val(1),
        mid, mv;
    do {
        mid = (l + r) / 2.0;
        mv = val(mid);
        if (Math.abs(mv) < 1e-6) {
            l = mid;
            lv = mv;
            break;
        } else if (sig(lv) === sig(mv)) {
            l = mid;
            lv = mv;
        } else {
            r = mid;
            rv = mv;
        }
    } while (rv > lv + 1e-6);
    var t = l;
    return (x - this.points[i].x) / (this.points[npi].x - this.points[i].x) * (this.variance[i + 1] - this.variance[i]) + this.variance[i] + (
        this._h00(t) * this._getvalue(i) +
        this._h10(t) * this.tangent[i][1] +
        this._h01(t) * this._getvalue(npi) +
        this._h11(t) * this.tangent[npi][1]);
};

CSpline.prototype.bezierPoints = function (i) {
    return [{
        x: this.points[i].x + this.tangent[i][0] / 3,
        y: this.points[i].y + this.tangent[i][1] / 3
    }, {
        x: this.points[i + 1].x - this.tangent[i + 1][0] / 3,
        y: this.points[i + 1].y - this.tangent[i + 1][1] / 3
    }];
};

CSpline.prototype.line = function (scaleX, scaleY, sampleNumber, low, high) {
    var scale = function (pt) {
            return {
                x: pt.x,
                y: scaleY(pt.y)
            };
        },
        cmd = pointString({
            x: scaleX(this.points[0].x),
            y: scaleY(this.points[0].y)
        });
    if (low === undefined) {
        low = scaleX(this.points[0].x + 1e-6);
    }
    low += 1e-6;
    if (high === undefined) {
        high = scaleX(lastOf(this.points).x - 1e-6);
    }
    high -= 1e-6;
    sampler(low, high, sampleNumber, function (x) {
        cmd += "L" + pointString(scale({x: x, y: this.evaluate(scaleX.invert(x))}));
    }, this);
    return cmd;
};

CSpline.prototype.lineReversed = function (scaleX, scaleY, sampleNumber, low, high) {
    var scale = function (pt) {
            return {
                x: pt.x,
                y: scaleY(pt.y)
            };
        },
        cmd = pointString({
            x: scaleX(lastOf(this.points).x),
            y: scaleY(lastOf(this.points).y)
        });
    if (low === undefined) {
        low = scaleX(this.points[0].x);
    }
    low += 1e-6;
    if (high === undefined) {
        high = scaleX(lastOf(this.points).x);
    }
    high -= 1e-6;
    sampler(high, low, sampleNumber, function (x) {
        cmd += "L" + pointString(scale({x: x, y: this.evaluate(scaleX.invert(x))}));
    }, this);
    return cmd;
};

CSpline.prototype.bezier = function (scaleX, scaleY) {
    var scale = function (pt) {
            return {
                x: scaleX(pt.x),
                y: scaleY(pt.y)
            };
        },
        cmd = pointString(scale(this.points[0])),
        prevX = this.points[0].x;
    for (var i = 1; i < this.points.length; i++) {
        if (this.points[i].x - prevX < 1e-6) {
            cmd += "L" + pointString(scale(this.points[i]));
        } else {
            var cp = this.bezierPoints(i - 1);
            cmd += "C" + pointString(scale(cp[0])) + "," +
                pointString(scale(cp[1])) + "," +
                pointString(scale(this.points[i]));
        }
        prevX = this.points[i].x;
    }
    return cmd;
};

CSpline.prototype.bezierReversed = function (scaleX, scaleY) {
    var scale = function (pt) {
            return {
                x: scaleX(pt.x),
                y: scaleY(pt.y)
            };
        },
        cmd = pointString(scale(lastOf(this.points))),
        prevX = lastOf(this.points).x;
    for (var i = this.points.length - 2; i >= 0; i--) {
        if (prevX - this.points[i].x < 1e-6) {
            cmd += "L" + pointString(scale(this.points[i]));
        } else {
            var cp = this.bezierPoints(i);
            cmd += "C" + pointString(scale(cp[1])) + "," +
                pointString(scale(cp[0])) + "," +
                pointString(scale(this.points[i]));
        }
        prevX = this.points[i].x;
    }
    return cmd;
};

function drawSpline(svg, xmin, xmax, f, scaleX, scaleY, n) {
    n = n || 200;
    var pts = [];
    sampler(xmin, xmax, n, function (x) {
        pts.push({
            x: x,
            y: f(x)
        });
    });
    drawPoints(svg, pts, scaleX, scaleY);
}

function drawSplineInBeizer(svg, spline, scaleX, scaleY) {
    svg.selectAll("_")
        .data([1])
        .enter().append("path")
        .attr("d", "M" + spline.bezier(scaleX, scaleY))
        .attr("stroke", "black")
        .attr("fill", "none")
        .attr("stroke-width", "1");
}

function testcspline(svg, data, scaleX, scaleY) {
    var layer = 1,
        pts = data[layer].map(function (d) {
            return {
                x: d.x,
                y: d.y0
            };
        });
    var spline = new CSpline(pts);
    //drawSpline(svg, firstOf(data[layer]).x, lastOf(data[layer]).x,
    //    function (x) { return spline.evaluate(x); }, scaleX, scaleY);
    drawSplineInBeizer(svg, spline, scaleX, scaleY);
}