Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /trunk/src/compiler/IL/kernel.sml
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log 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