(* 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
(* kernel name *)
val name : kernel -> string
(* kernel support *)
val support : kernel -> int
(* representation of i'th derivative of the kernel *)
val curve : kernel * int -> {
isOdd : bool,
isCont : bool,
segs : polynomial list (* piece-wise polynomial that defines *)
(* the curve over the positive support *)
}
val evaluate : polynomial * int -> Rational.rat
(* some standard kernels *)
val tent : kernel (* linear interpolation *)
val ctmr : kernel (* Catmull-Rom interpolation *)
val bspl3 : kernel (* cubic bspline reconstruction, doesn't interpolate *)
val bspl5 : kernel (* quintic bspline reconstruction, doesn't interpolate *)
end = struct
structure R = Rational
structure A = Array
val maxDiffLevels = 15 (* support upto 15 levels of differentiation *)
type coefficient = R.rat
val zero = R.zero
val one = R.fromInt 1
(* polynomial represented as list of coefficients, where ith element is
* coefficient for x^i.
*)
type polynomial = coefficient list
fun differentiate [] = raise Fail "invalid polynomial"
| differentiate [_] = [zero]
| differentiate (_::coeffs) = let
fun lp (_, []) = []
| lp (i, c::r) = R.*(R.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 = R.fromInt x
fun eval (sum, [], xn) = sum
| eval (sum, c::r, xn) = eval(R.+(sum, R.*(c, xn)), r, R.*(x, xn))
in
eval (zero, poly, one)
end
type curve = {
isOdd : bool,
isCont : bool,
segs : polynomial list (* piece-wise polynomial that defines *)
(* the curve over the positive support *)
}
datatype kernel = K of {
name : string,
support : int, (* number of samples to left/right *)
curves : curve option array (* cache of curves indexed by differentiation level *)
}
(* determine if a list of polynomials represents a continuous piece-wise polynomial *)
fun isContinuous [_] = true
| isContinuous (f0::r) = let
fun chk (i, f_i, []) = (R.zero = evaluate(f_i, i))
| chk (i, f_i, f_i1::r) = let
val y_i = evaluate(f_i, i)
val y_i1 = evaluate(f_i1, i)
in
if (y_i = y_i1)
then chk(i+1, f_i1, r)
else
end
in
chk (1, f0, r)
end
(* kernel name *)
fun name (K{name, ...}) = name
(* kernel support *)
fun support (K{support, ...}) = support
(* representation of i'th derivative of the kernel *)
fun curve (K{curves, ...}, k) = (case A.sub(curves, k)
of SOME curve => curve
| NONE => let
(* compute the (k+1)'th derivative, given the k'th *)
fun diff (k, {isOdd, isCont, segs}) = let
val segs' = List.map differentiate segs
val isOdd' = not isOdd
val isCont' = if isCont andalso isContinuous segs'
then (not isOdd') orelse (evaluate(List.hd segs', 0) = R.zero)
else false
in {
isOdd = isOdd',
isCont = isCont',
segs = segs'
} end
fun lp (j, curve) = if (j < k)
then (case A.sub(curves, j+1)
of NONE => let
val curve' = diff(j+1, curve)
in
A.update(curves, j+1, SOME curve');
lp (j+1, curve')
end
| SOME curve' => lp(j+1, curve')
(* end case *))
else curve
in
lp (0, valOf(A.sub(curves, 0)))
end
(* end case *))
(* some standard kernels *)
local
val op / = R./
fun r i = R.fromInt i
fun mkKernel {name, support, segs} = let
val curves = Array.array(maxDiffLevels+1, NONE)
val curve0 = {
isOdd = false,
isCont = isContinuous segs,
segs = segs
}
in
A.update (curves, 0, SOME curve0);
K{name=name, support=support, curves=curves}
end
in
val tent : kernel = mkKernel{ (* linear interpolation *)
name = "tent",
support = 1,
segs = [[r 1, r ~1]]
}
val ctmr : kernel = mkKernel{ (* Catmull-Rom interpolation *)
name = "ctmr",
support = 2,
segs = [
[r 1, r 0, ~5/2, 3/2],
[r 2, r ~4, 5/2, ~1/2]
]
}
val bspl3 : kernel = mkKernel{ (* cubic bspline reconstruction, doesn't interpolate *)
name = "bspl3",
support = 2,
segs = [
[ 2/3, r 0, r ~1, 1/2 ],
[ 4/3, r ~2, r 1, ~1/6 ]
]
}
val bspl5 : kernel = mkKernel{ (* quintic bspline reconstruction, doesn't interpolate *)
name = "bspl5",
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