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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 435 - (view) (download)

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

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