Cardinal splines in Scala

splines.scala

class Point(val x: Double, val y : Double) {
}

object Point {
  def apply(x : Double, y : Double) = new Point(x, y)
  def unapply(p : Point) = Some((p.x, p.y))
}

object CardinalSpline {
  import scalala.tensor.dense._

  // tension factor; at 0 we have Catmull-Rom splines
  val t = 0.0
  
  // The number of straight lines used to approximate each segment of
  // the curve
  val steps = 20

  val s = (1-t)/2

  val basis = DenseMatrix((0.0, 1.0, 0.0,   0.0),
                          ( -s, 0.0, s,     0.0),
                          (2*s, s-3, 3-2*s, -s),
                          (-s,  2-s, s-2,   s))

  // precompute the transformation matrices
  val uBs = (0.0 to 1.0 by (1.0/steps)).view.map {u =>
      val U = DenseVector(1, u, u*u, u*u*u)
      U.t * basis
  }
  
  def segment(ps : Array[Point], i : Int) = {
    val px = DenseVector(ps(i).x, ps(i+1).x, ps(i+2).x, ps(i+3).x)
    val py = DenseVector(ps(i).y, ps(i+1).y, ps(i+2).y, ps(i+3).y)
    uBs.map(uB => Point(px dot uB, py dot uB))
  }

  def path(ps : Array[Point]) = {
    // we compute the endpoints by simple reflection
    val n = ps.size
    val p0 = Point(2 * ps(0).x - ps(1).x, 2 * ps(0).y - ps(1).y)
    val p1 = Point(2 * ps(n-1).x - ps(n-2).x,
                   2 * ps(n-1).y - ps(n-2).y)
    // and compute each segment piecewise
    val points = p0 +: ps :+ p1
    (0 until (n - 1)).flatMap(segment(points, _))
  }
}