(* kernel.sml
 *
 * COPYRIGHT (c) 2010 The Diderot Project (http://diderot.cs.uchicago.edu)
 * All rights reserved.
 *
 * QUESTION: should we 
 *)

structure Kernel : sig

    type coefficient = Rational.rat

  (* polynomial represented as list of coefficients, where ith element is
   * coefficient for x^i.
   *)
    type polynomial = coefficient list

    type kernel = {
	support : int,		(* number of samples to left/right *)
	segs : polynomial list,	(* piece-wise polynomial that defines the kernel.  Since *)
				(* a kernel is symmetric, we only store the positive *)
				(* segments, with segs[i] covering the domain [i,i+1) *)
      }

    val differentiate : polynomial -> polynomial

    val evaluate : polynomial * int -> Rational.rat

  end = struct

    structure R = Rational

    type coefficient = R.rat

    val zero = R.fromInt 0
    val one = R.fromInt 1

  (* polynomial represented as list of coefficients, where ith element is
   * coefficient for x^i.
   *)
    type polynomial = coefficient list

    type kernel = {
	support : int,		(* number of samples to left/right *)
	segs : polynomial list,	(* piece-wise polynomial that defines the *)
				(* kernel; there are 2*support segments *)
      }

    fun differentiate [] = raise Fail "invalid polynomial"
      | differentiate [_] = [zero]
      | differentiate (_::coeffs) = let
	  fun lp (_, []) = []
	    | lp (i, c::r) = R.mul(F.fromInt i, c) :: lp(i+1, r)
	  in
	    lp (1, coeffs)
	  end

  (* evaluate a polynomial at an integer coordinate (used to test continuity) *)
    fun evaluate (poly, x) = let
	  val x = F.fromInt x
	  fun eval (sum, [], xn) = sum
	    | eval (sum, c::r, xn) = eval(F.add(sum, R.mul(c, xn)), r, R.mul(x, xn))
	  in
	    eval (zero, poly, one)
	  end

  (* some standard kernels *)
    local
      val op / = R./
      fun r i = R.fromInt i
    in
    val tent : kernel = {	(* linear interpolation *)
	    support = 1,
	    segs = [[r 1, r ~1]]
	  }
    val ctmr : kernel = {	(* Catmull-Rom interpolation *)
	    support = 2,
	    segs = [
		[r 1, r 0, ~5/2, 3/2],
		[r 2, r ~4, 5/2, ~1/2]
	      ]
	  }
    val bspl3 : kernel = {	(* cubic bspline reconstruction, doesn't interpolate *)
	    support = 2,
	    segs = [
		[ 2/3, r 0, r ~1, 1/2 ],
		[ 4/3, r ~2, r 1, ~1/6 ]'
	      ]
	  }
    val bspl5 : kernel = {	(* quintic bspline reconstruction, doesn't interpolate *)
	    support = 3,
	    segs = [
		[ 11/20, r 0, ~1/2, r 0, 1/4, ~1/12 ],
		[ 17/40, 5/8, ~7/4, 5/4, ~3/8, 1/24 ],
		[ 81/40, ~27/8, 9/4, ~3/4, 1/8, ~1/120 ]
	      ]
	  }
    end

  end