/* calculates the natural cubic spline that interpolates
y[0], y[1], ... y[n]
The first segment is returned as
C[0].a + C[0].b*u + C[0].c*u^2 + C[0].d*u^3 0<=u <1
the other segments are in C[1], C[2], ...  C[n-1]
438146
public Cubic[] calcNaturalCubic( int n, double[] x ) {
        double[] gamma = new double[n + 1];
        double[] delta = new double[n + 1];
        double[] D = new double[n + 1];
        int i;
        /* We solve the equation
           [2 1       ] [D[0]]   [3(x[1] - x[0])  ]
           |1 4 1     | |D[1]|   |3(x[2] - x[0])  |
           |  1 4 1   | | .  | = |      .         |
           |    ..... | | .  |   |      .         |
           |     1 4 1| | .  |   |3(x[n] - x[n-2])|
           [       1 2] [D[n]]   [3(x[n] - x[n-1])]
           
           by using row operations to convert the matrix to upper triangular
           and then back sustitution.  The D[i] are the derivatives at the knots.
           */

        gamma[0] = 1.0f / 2.0f;
        for( i = 1; i < n; i++ ) {
            gamma[i] = 1 / (4 - gamma[i - 1]);
        }
        gamma[n] = 1 / (2 - gamma[n - 1]);

        delta[0] = 3 * (x[1] - x[0]) * gamma[0];
        for( i = 1; i < n; i++ ) {
            delta[i] = (3 * (x[i + 1] - x[i - 1]) - delta[i - 1]) * gamma[i];
        }
        delta[n] = (3 * (x[n] - x[n - 1]) - delta[n - 1]) * gamma[n];

        D[n] = delta[n];
        for( i = n - 1; i >= 0; i-- ) {
            D[i] = delta[i] - gamma[i] * D[i + 1];
        }

        /* now compute the coefficients of the cubics */
        Cubic[] C = new Cubic[n];
        for( i = 0; i < n; i++ ) {
            C[i] = new Cubic((double) x[i], D[i], 3 * (x[i + 1] - x[i]) - 2 * D[i] - D[i + 1], 2 * (x[i] - x[i + 1]) + D[i]
                    + D[i + 1]);
        }
        return C;
    }
149480
def interpolation_linear(x, x1, x2, y1, y2):
    
    m = (y2 - y1) / (x2 - x1)
    t = (x - x1)
    return m * t + y1
215808
function (x1, y1, x2, y2, x3, y3) {
        var dashSum = this._dashSum;
        var offset = this._dashOffset;
        var lineDash = this._lineDash;
        var ctx = this._ctx;

        var x0 = this._xi;
        var y0 = this._yi;
        var t;
        var dx;
        var dy;
        var cubicAt = curve.cubicAt;
        var bezierLen = 0;
        var idx = this._dashIdx;
        var nDash = lineDash.length;

        var x;
        var y;

        var tmpLen = 0;

        if (offset < 0) {
            // Convert to positive offset
            offset = dashSum + offset;
        }
        offset %= dashSum;
        // Bezier approx length
        for (t = 0; t < 1; t += 0.1) {
            dx = cubicAt(x0, x1, x2, x3, t + 0.1)
                - cubicAt(x0, x1, x2, x3, t);
            dy = cubicAt(y0, y1, y2, y3, t + 0.1)
                - cubicAt(y0, y1, y2, y3, t);
            bezierLen += mathSqrt(dx * dx + dy * dy);
        }

        // Find idx after add offset
        for (; idx < nDash; idx++) {
            tmpLen += lineDash[idx];
            if (tmpLen > offset) {
                break;
            }
        }
        t = (tmpLen - offset) / bezierLen;

        while (t <= 1) {

            x = cubicAt(x0, x1, x2, x3, t);
            y = cubicAt(y0, y1, y2, y3, t);

            // Use line to approximate dashed bezier
            // Bad result if dash is long
            idx % 2 ? ctx.moveTo(x, y)
                : ctx.lineTo(x, y);

            t += lineDash[idx] / bezierLen;

            idx = (idx + 1) % nDash;
        }

        // Finish the last segment and calculate the new offset
        (idx % 2 !== 0) && ctx.lineTo(x3, y3);
        dx = x3 - x;
        dy = y3 - y;
        this._dashOffset = -mathSqrt(dx * dx + dy * dy);
    }
264365
def linear_interpolate(x1, y1, x2):
    
    li = LinInt(x1, y1)
    return li(x2)
737469
function (bounds, vertex, next, orthogonal)
	{
		var x = bounds.x;
		var y = bounds.y;
		var w = bounds.width;
		var h = bounds.height;
		
		var cx = x + w / 2;
		var cy = y + h / 2;

		var px = next.x;
		var py = next.y;

		// Special case for intersecting the diamond's corners
		if (cx == px)
		{
			if (cy > py)
			{
				return new mxPoint(cx, y); // top
			}
			else
			{
				return new mxPoint(cx, y + h); // bottom
			}
		}
		else if (cy == py)
		{
			if (cx > px)
			{
				return new mxPoint(x, cy); // left
			}
			else
			{
				return new mxPoint(x + w, cy); // right
			}
		}
		
		var tx = cx;
		var ty = cy;
		
		if (orthogonal)
		{
			if (px >= x && px <= x + w)
			{
				tx = px;
			}
			else if (py >= y && py <= y + h)
			{
				ty = py;
			}
		}
		
		// In which quadrant will the intersection be?
		// set the slope and offset of the border line accordingly
		if (px < cx)
		{
			if (py < cy)
			{
				return mxUtils.intersection(px, py, tx, ty, cx, y, x, cy);
			}
			else
			{
				return mxUtils.intersection(px, py, tx, ty, cx, y + h, x, cy);
			}
		}
		else if (py < cy)
		{
			return mxUtils.intersection(px, py, tx, ty, cx, y, x + w, cy);
		}
		else
		{
			return mxUtils.intersection(px, py, tx, ty, cx, y + h, x + w, cy);
		}
	}
576421
function (x0, y0, x1, y1, color, c)
    {
        PixelArt.put(x, y, color, c)
        var dx = x1 - x0;
        var dy = y1 - y0;
        var adx = Math.abs(dx);
        var ady = Math.abs(dy);
        var eps = 0;
        var sx = dx > 0 ? 1 : -1;
        var sy = dy > 0 ? 1 : -1;
        if (adx > ady)
        {
            for (var x = x0, y = y0; sx < 0 ? x >= x1 : x <= x1; x += sx)
            {
                PixelArt.put(x, y)
                eps += ady;
                if ((eps << 1) >= adx)
                {
                    y += sy;
                    eps -= adx;
                }
            }
        } else
        {
            for (var x = x0, y = y0; sy < 0 ? y >= y1 : y <= y1; y += sy)
            {
                PixelArt.put(x, y)
                eps += adx;
                if ((eps << 1) >= ady)
                {
                    x += sx;
                    eps -= ady;
                }
            }
        }
    }
890878
def to_curve_spline(obj):
    '''
    
    '''
    if   is_curve_spline(obj):            return obj
    elif is_tuple(obj) and len(obj) == 2: (crds,opts) = obj
    else:                                 (crds,opts) = (obj,{})
    if pimms.is_matrix(crds) or is_curve_spline(crds): crds = [crds]
    spls = [c for c in crds if is_curve_spline(c)]
    opts = dict(opts)
    if 'weights' not in opts and len(spls) == len(crds):
        if all(c.weights is not None for c in crds):
            opts['weights'] = np.concatenate([c.weights for c in crds])
    if 'order' not in opts and len(spls) > 0:
        opts['order'] = np.min([c.order for c in spls])
    if 'smoothing' not in opts and len(spls) > 0:
        sm = set([c.smoothing for c in spls])
        if len(sm) == 1: opts['smoothing'] = list(sm)[0]
        else: opts['smoothing'] = None
    crds = [x.crds if is_curve_spline(crds) else np.asarray(x) for x in crds]
    crds = [x if x.shape[0] == 2 else x.T for x in crds]
    crds = np.hstack(crds)
    return curve_spline(crds, **opts)
587718
def get_Zstar(x_segment, y_segment, x_int, y_int, slope, n):
    
    total = 0
    first_time = True
    for num, x in enumerate(x_segment):
        b_wiggle = get_b_wiggle(x, y_segment[num], y_int)
        result = 100 * ( old_div((x * abs(b_wiggle - abs(slope)) ), abs(y_int)) )
        total += result
        first_time = False
    Zstar = (old_div(1., (n - 1.))) * total
    return Zstar
1341556
public int[] makeLut() {
        int numKnots = x.length;
        float[] nx = new float[numKnots + 2];
        float[] ny = new float[numKnots + 2];
        System.arraycopy(x, 0, nx, 1, numKnots);
        System.arraycopy(y, 0, ny, 1, numKnots);
        nx[0] = nx[1];
        ny[0] = ny[1];
        nx[numKnots + 1] = nx[numKnots];
        ny[numKnots + 1] = ny[numKnots];

        int[] table = new int[256];
        for (int i = 0; i < 1024; i++) {
            float f = i / 1024.0f;
            int x = (int) (255 * Curve.Spline(f, nx.length, nx) + 0.5f);
            int y = (int) (255 * Curve.Spline(f, nx.length, ny) + 0.5f);
            x = x > 255 ? 255 : x;
            x = x < 0 ? 0 : x;
            y = y > 255 ? 255 : y;
            y = y < 0 ? 0 : y;
            table[x] = y;
        }
        return table;
    }
820467
def get_Z(x_segment, y_segment, x_int, y_int, slope):
    
    Z = 0
    first_time = True
    for num, x in enumerate(x_segment):
        b_wiggle = get_b_wiggle(x, y_segment[num], y_int)
        z = old_div((x * abs(b_wiggle - abs(slope)) ), abs(x_int))
        Z += z
        first_time = False
    return Z
505101
def interpolate_cubic(self, lons, lats, data):
        
        return self.interpolate(lons, lats, data, order=3)
54131
function () {
            var i;
            var edges = this.edges;
            var normals = this.normals;
            var indices = this.indices;

            // Copy the original points array and apply the offset/angle
            var points = this.points;
            var len = points.length;

            if (len < 3) {
                throw new Error("Requires at least 3 points");
            }

            // Calculate the edges/normals
            for (i = 0; i < len; i++) {
                if (edges[i] === undefined) {
                    edges[i] = new me.Vector2d();
                }
                edges[i].copy(points[(i + 1) % len]).sub(points[i]);

                if (normals[i] === undefined) {
                    normals[i] = new me.Vector2d();
                }
                normals[i].copy(edges[i]).perp().normalize();
            }
            // trunc array
            edges.length = len;
            normals.length = len;
            // do not do anything here, indices will be computed by
            // toIndices if array is empty upon function call
            indices.length = 0;

            return this;
        }
525257
private static double interpolateUsingXArrAndYStride(final double[] xArr, final double yStride,
      final int offset, final double x) {
    return cubicInterpolate(
        xArr[offset + 0], yStride * (offset + 0),
        xArr[offset + 1], yStride * (offset + 1),
        xArr[offset + 2], yStride * (offset + 2),
        xArr[offset + 3], yStride * (offset + 3),
        x);
  }
294778
function curve($x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3) {
    // in the current line style, draw a bezier curve from (x0,y0) to (x3,y3) using the other two points
    // as the control points for the curve.
    $this->objects[$this->currentContents]['c'] .=
      sprintf("\n%.3F %.3F m %.3F %.3F %.3F %.3F %.3F %.3F c S", $x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3);
  }
516862
function(x0, y0, x1, y1, x2, y2, x3, y3) {
  var tvalues = [],
      bounds = [[], []],
      a, b, c, t, t1, t2, b2ac, sqrtb2ac;
  for (var i = 0; i < 2; ++i) {
    if (i == 0) {
      b = 6 * x0 - 12 * x1 + 6 * x2;
      a = -3 * x0 + 9 * x1 - 9 * x2 + 3 * x3;
      c = 3 * x1 - 3 * x0;
    } else {
      b = 6 * y0 - 12 * y1 + 6 * y2;
      a = -3 * y0 + 9 * y1 - 9 * y2 + 3 * y3;
      c = 3 * y1 - 3 * y0;
    }
    if (Math.abs(a) < 1e-12) {
      if (Math.abs(b) < 1e-12) {
        continue;
      }
      t = -c / b;
      if (0 < t && t < 1) {
        tvalues.push(t);
      }
      continue;
    }
    b2ac = b * b - 4 * c * a;
    sqrtb2ac = Math.sqrt(b2ac);
    if (b2ac < 0) {
      continue;
    }
    t1 = (-b + sqrtb2ac) / (2 * a);
    if (0 < t1 && t1 < 1) {
      tvalues.push(t1);
    }
    t2 = (-b - sqrtb2ac) / (2 * a);
    if (0 < t2 && t2 < 1) {
      tvalues.push(t2);
    }
  }

  var x, y, j = tvalues.length,
      jlen = j,
      mt;
  while (j--) {
    t = tvalues[j];
    mt = 1 - t;
    bounds[0][j] = (mt * mt * mt * x0) + (3 * mt * mt * t * x1) + (3 * mt * t * t * x2) + (t * t * t * x3);
    bounds[1][j] = (mt * mt * mt * y0) + (3 * mt * mt * t * y1) + (3 * mt * t * t * y2) + (t * t * t * y3);
  }

  bounds[0][jlen] = x0;
  bounds[1][jlen] = y0;
  bounds[0][jlen + 1] = x3;
  bounds[1][jlen + 1] = y3;
  bounds[0].length = bounds[1].length = jlen + 2;

  return {
    min: {x: Math.min.apply(0, bounds[0]), y: Math.min.apply(0, bounds[1])},
    max: {x: Math.max.apply(0, bounds[0]), y: Math.max.apply(0, bounds[1])}
  };
}
1008607
function(d2) {

      var w = this['w'];
      var x = this['x'];
      var y = this['y'];
      var z = this['z'];

      var n = w * w + x * x + y * y + z * z;
      var s = n === 0 ? 0 : 2 / n;
      var wx = s * w * x, wy = s * w * y, wz = s * w * z;
      var xx = s * x * x, xy = s * x * y, xz = s * x * z;
      var yy = s * y * y, yz = s * y * z, zz = s * z * z;

      if (d2) {
        return [
          [1 - (yy + zz), xy - wz, xz + wy],
          [xy + wz, 1 - (xx + zz), yz - wx],
          [xz - wy, yz + wx, 1 - (xx + yy)]];
      }

      return [
        1 - (yy + zz), xy - wz, xz + wy,
        xy + wz, 1 - (xx + zz), yz - wx,
        xz - wy, yz + wx, 1 - (xx + yy)];
    }
854983
def interpolation_points(self, N):
        
        if N == 1:
            return np.array([0.])
        return np.cos(np.arange(N)*np.pi/(N-1))