Annotation of /trunk/src/compiler/IL/kernel.sml

Parent Directory | Revision Log

---

Revision 155 - (view) (download)

1 :	jhr	108	(* kernel.sml
2 :			*
3 :			* COPYRIGHT (c) 2010 The Diderot Project (http://diderot.cs.uchicago.edu)
4 :			* All rights reserved.
5 :	jhr	117	*
6 :			* QUESTION: should we
7 :	jhr	108	*)
8 :
9 :			structure Kernel : sig
10 :
11 :	jhr	117	type coefficient = Rational.rat
12 :	jhr	108
13 :			(* polynomial represented as list of coefficients, where ith element is
14 :			* coefficient for x^i.
15 :			*)
16 :			type polynomial = coefficient list
17 :
18 :	jhr	139	type kernel
19 :	jhr	108
20 :	jhr	139	(* kernel name *)
21 :			val name : kernel -> string
22 :	jhr	108
23 :	jhr	139	(* kernel support *)
24 :			val support : kernel -> int
25 :
26 :			(* representation of i'th derivative of the kernel *)
27 :			val curve : kernel * int -> {
28 :			isOdd : bool,
29 :			isCont : bool,
30 :			segs : polynomial list (* piece-wise polynomial that defines *)
31 :			(* the curve over the positive support *)
32 :			}
33 :
34 :	jhr	117	val evaluate : polynomial * int -> Rational.rat
35 :	jhr	108
36 :	jhr	153	(* some standard kernels *)
37 :			val tent : kernel (* linear interpolation *)
38 :			val ctmr : kernel (* Catmull-Rom interpolation *)
39 :			val bspl3 : kernel (* cubic bspline reconstruction, doesn't interpolate *)
40 :			val bspl5 : kernel (* quintic bspline reconstruction, doesn't interpolate *)
41 :
42 :	jhr	108	end = struct
43 :
44 :	jhr	117	structure R = Rational
45 :	jhr	139	structure A = Array
46 :	jhr	108
47 :	jhr	139	val maxDiffLevels = 15 (* support upto 15 levels of differentiation *)
48 :
49 :	jhr	117	type coefficient = R.rat
50 :	jhr	108
51 :	jhr	153	val zero = R.zero
52 :	jhr	117	val one = R.fromInt 1
53 :
54 :	jhr	108	(* polynomial represented as list of coefficients, where ith element is
55 :			* coefficient for x^i.
56 :			*)
57 :			type polynomial = coefficient list
58 :
59 :			fun differentiate [] = raise Fail "invalid polynomial"
60 :	jhr	117	| differentiate [_] = [zero]
61 :	jhr	108	| differentiate (_::coeffs) = let
62 :			fun lp (_, []) = []
63 :	jhr	139	| lp (i, c::r) = R.*(R.fromInt i, c) :: lp(i+1, r)
64 :	jhr	108	in
65 :			lp (1, coeffs)
66 :			end
67 :
68 :			(* evaluate a polynomial at an integer coordinate (used to test continuity) *)
69 :			fun evaluate (poly, x) = let
70 :	jhr	139	val x = R.fromInt x
71 :	jhr	108	fun eval (sum, [], xn) = sum
72 :	jhr	139	| eval (sum, c::r, xn) = eval(R.+(sum, R.*(c, xn)), r, R.*(x, xn))
73 :	jhr	108	in
74 :	jhr	117	eval (zero, poly, one)
75 :	jhr	108	end
76 :
77 :	jhr	139	type curve = {
78 :			isOdd : bool,
79 :			isCont : bool,
80 :			segs : polynomial list (* piece-wise polynomial that defines *)
81 :			(* the curve over the positive support *)
82 :			}
83 :
84 :			datatype kernel = K of {
85 :			name : string,
86 :			support : int, (* number of samples to left/right *)
87 :			curves : curve option array (* cache of curves indexed by differentiation level *)
88 :			}
89 :
90 :			(* determine if a list of polynomials represents a continuous piece-wise polynomial *)
91 :	jhr	153	fun isContinuous [_] = true
92 :			| isContinuous (f0::r) = let
93 :	jhr	139	fun chk (i, f_i, []) = (R.zero = evaluate(f_i, i))
94 :			| chk (i, f_i, f_i1::r) = let
95 :			val y_i = evaluate(f_i, i)
96 :			val y_i1 = evaluate(f_i1, i)
97 :			in
98 :	jhr	155	(y_i = y_i1) andalso chk(i+1, f_i1, r)
99 :	jhr	139	end
100 :			in
101 :	jhr	153	chk (1, f0, r)
102 :	jhr	139	end
103 :
104 :			(* kernel name *)
105 :			fun name (K{name, ...}) = name
106 :
107 :			(* kernel support *)
108 :			fun support (K{support, ...}) = support
109 :
110 :			(* representation of i'th derivative of the kernel *)
111 :			fun curve (K{curves, ...}, k) = (case A.sub(curves, k)
112 :			of SOME curve => curve
113 :			| NONE => let
114 :			(* compute the (k+1)'th derivative, given the k'th *)
115 :			fun diff (k, {isOdd, isCont, segs}) = let
116 :			val segs' = List.map differentiate segs
117 :	jhr	153	val isOdd' = not isOdd
118 :			val isCont' = if isCont andalso isContinuous segs'
119 :			then (not isOdd') orelse (evaluate(List.hd segs', 0) = R.zero)
120 :			else false
121 :	jhr	139	in {
122 :	jhr	153	isOdd = isOdd',
123 :			isCont = isCont',
124 :	jhr	139	segs = segs'
125 :			} end
126 :			fun lp (j, curve) = if (j < k)
127 :			then (case A.sub(curves, j+1)
128 :			of NONE => let
129 :			val curve' = diff(j+1, curve)
130 :			in
131 :			A.update(curves, j+1, SOME curve');
132 :			lp (j+1, curve')
133 :			end
134 :			| SOME curve' => lp(j+1, curve')
135 :			(* end case *))
136 :			else curve
137 :			in
138 :			lp (0, valOf(A.sub(curves, 0)))
139 :			end
140 :			(* end case *))
141 :
142 :	jhr	117	(* some standard kernels *)
143 :			local
144 :			val op / = R./
145 :			fun r i = R.fromInt i
146 :	jhr	139	fun mkKernel {name, support, segs} = let
147 :			val curves = Array.array(maxDiffLevels+1, NONE)
148 :			val curve0 = {
149 :			isOdd = false,
150 :			isCont = isContinuous segs,
151 :			segs = segs
152 :			}
153 :			in
154 :			A.update (curves, 0, SOME curve0);
155 :			K{name=name, support=support, curves=curves}
156 :			end
157 :	jhr	117	in
158 :	jhr	139	val tent : kernel = mkKernel{ (* linear interpolation *)
159 :	jhr	135	name = "tent",
160 :	jhr	117	support = 1,
161 :			segs = [[r 1, r ~1]]
162 :			}
163 :	jhr	139	val ctmr : kernel = mkKernel{ (* Catmull-Rom interpolation *)
164 :	jhr	135	name = "ctmr",
165 :	jhr	117	support = 2,
166 :			segs = [
167 :			[r 1, r 0, ~5/2, 3/2],
168 :			[r 2, r ~4, 5/2, ~1/2]
169 :			]
170 :			}
171 :	jhr	139	val bspl3 : kernel = mkKernel{ (* cubic bspline reconstruction, doesn't interpolate *)
172 :	jhr	135	name = "bspl3",
173 :	jhr	117	support = 2,
174 :			segs = [
175 :	jhr	122	[ 2/3, r 0, r ~1, 1/2 ],
176 :	jhr	135	[ 4/3, r ~2, r 1, ~1/6 ]
177 :	jhr	117	]
178 :			}
179 :	jhr	139	val bspl5 : kernel = mkKernel{ (* quintic bspline reconstruction, doesn't interpolate *)
180 :	jhr	135	name = "bspl5",
181 :	jhr	117	support = 3,
182 :			segs = [
183 :			[ 11/20, r 0, ~1/2, r 0, 1/4, ~1/12 ],
184 :	jhr	127	[ 17/40, 5/8, ~7/4, 5/4, ~3/8, 1/24 ],
185 :	jhr	117	[ 81/40, ~27/8, 9/4, ~3/4, 1/8, ~1/120 ]
186 :			]
187 :			}
188 :			end
189 :
190 :	jhr	108	end

---

root@smlnj-gforge.cs.uchicago.edu	ViewVC Help
Powered by ViewVC 1.0.0