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

SCM Repository

[diderot] Annotation of /benchmarks/programs/lic2d/bmark-teem.c
ViewVC logotype

Annotation of /benchmarks/programs/lic2d/bmark-teem.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1612 - (view) (download) (as text)

1 : nseltzer 1592 /*! \file /lic2d/bmark-teem.c
2 :     *
3 : glk 1612 * \author Nick Seltzer and Gordon Kindlmann
4 : nseltzer 1592 */
5 :    
6 :     /*
7 :     * COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
8 :     * All rights reserved.
9 :     */
10 :    
11 :     #include <stdio.h>
12 :     #include <math.h>
13 :     #include <stdbool.h>
14 :     #include "teem/nrrd.h"
15 : glk 1612 #include "teem/ell.h"
16 : nseltzer 1592 #include "teem/gage.h"
17 :    
18 : jhr 1596 #include "teem-defs.h"
19 : nseltzer 1592
20 :     STATIC_INLINE void scale2(double scl, double vec[2], double res[2])
21 :     {
22 :     res[0] = scl * vec[0];
23 :     res[1] = scl * vec[1];
24 :     }
25 :    
26 :     STATIC_INLINE void add2(double vec1[2], double vec2[2], double res[2])
27 :     {
28 :     res[0] = vec1[0] + vec2[0];
29 :     res[1] = vec1[1] + vec2[1];
30 :     }
31 :    
32 :     STATIC_INLINE double dot2(double vec1[2], double vec2[2])
33 :     {
34 :     return (vec1[0] * vec2[0] + vec1[1] * vec2[1]);
35 :     }
36 :    
37 :     STATIC_INLINE double mag2(double vec[2])
38 :     {
39 :     return sqrt(dot2(vec, vec));
40 :     }
41 :    
42 :     STATIC_INLINE void norm2(double vec[2], double res[2])
43 :     {
44 :     double mag = mag2(vec);
45 :     res[0] = vec[0] / mag;
46 :     res[1] = vec[1] / mag;
47 :     }
48 :    
49 :     STATIC_INLINE bool inside(double pos[2])
50 :     {
51 :     // XXX - Hack
52 : nseltzer 1608 if (pos[0] < 0.04 || pos[0] > 6.76) return false;
53 :     if (pos[1] < 0.04 || pos[1] > 3.7) return false;
54 : nseltzer 1592 return true;
55 :     }
56 :    
57 :     void atR(void *_data, double pos[2], double *out)
58 :     {
59 :     float *data = (float*)_data;
60 : nseltzer 1608 double temp[2][2];
61 : nseltzer 1592 int a = floor(pos[0] * 150);
62 :     int b = floor(pos[1] * 150);
63 : nseltzer 1608 for (int i = 0; i < 4; i++) {
64 :     for (int j = 0; j < 2; j++) {
65 :     temp[i][j] = *(data + (i + a) + (j + b) * 1020);
66 :     }
67 :     }
68 :     double temp2[2];
69 :     double t = pos[0]*150 - floor(pos[0]*150);
70 :     for (int i = 0; i < 2; i++) {
71 :     temp2[i] = (temp[0][i] * (1 - t) + temp[1][i] * t);
72 :     }
73 :     t = pos[1]*150 - floor(pos[1]*150);
74 :     *out = (temp2[0] * (1 - t) + temp2[1] * t);
75 : nseltzer 1592 }
76 :    
77 :     void at(void* _data, double pos[2], double out[2])
78 :     {
79 :     // XXX - Hack
80 :     float *data = (float*)_data;
81 :     double tempx[4][4];
82 :     double tempy[4][4];
83 :     int a = floor(pos[0] * 50) + 5;
84 :     int b = floor(pos[1] * 50) + 5;
85 : nseltzer 1608 for (int i = 0; i < 4; i++) {
86 :     for (int j = 0; j < 4; j++) {
87 : nseltzer 1592 tempx[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 0);
88 :     tempy[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 1);
89 :     }
90 :     }
91 : nseltzer 1608 double t = pos[0]*150 - floor(pos[0]*150);
92 : nseltzer 1592 double tempxx[4];
93 :     double tempyx[4];
94 : nseltzer 1608 for (int i = 0; i < 4; i++) {
95 : nseltzer 1592 tempxx[i] = 0.5 * (2 * tempx[0][i] +
96 :     (-tempx[0][i] + tempx[2][i]) * t +
97 :     (2 * tempx[0][i] - 5 * tempx[1][i] + 4 * tempx[2][i] - tempx[3][i]) * t * t +
98 :     (-tempx[0][i] + 3 * tempx[1][i] - 3 * tempx[2][i] + tempx[3][i]) * t * t * t);
99 :     tempyx[i] = 0.5 * (2 * tempy[0][i] +
100 :     (-tempy[0][i] + tempy[2][i]) * t +
101 :     (2 * tempy[0][i] - 5 * tempy[1][i] + 4 * tempy[2][i] - tempy[3][i]) * t * t +
102 :     (-tempy[0][i] + 3 * tempy[1][i] - 3 * tempy[2][i] + tempy[3][i]) * t * t * t);
103 :     }
104 : nseltzer 1608 t = pos[1]*150 - floor(pos[1]*150);
105 : nseltzer 1592 out[0] = 0.5 * (2 * tempxx[0] +
106 :     (-tempxx[0] + tempxx[2]) * t +
107 :     (2 * tempxx[0] - 5 * tempxx[1] + 4 * tempxx[2] - tempxx[3]) * t * t +
108 :     (-tempxx[0] + 3 * tempxx[1] - 3 * tempxx[2] + tempxx[3]) * t * t * t);
109 :     out[1] = 0.5 * (2 * tempyx[0] +
110 :     (-tempyx[0] + tempyx[2]) * t +
111 :     (2 * tempyx[0] - 5 * tempyx[1] + 4 * tempyx[2] - tempyx[3]) * t * t +
112 :     (-tempyx[0] + 3 * tempyx[1] - 3 * tempyx[2] + tempyx[3]) * t * t * t);
113 :     return;
114 :     }
115 :    
116 :     double crl_at(void* _data, double pos[2])
117 :     {
118 :     // XXX - Hack
119 :     float *data = (float*)_data;
120 :     double tempx[4][4];
121 :     double tempy[4][4];
122 :     double temp[2];
123 :     int a = floor(pos[0] * 50) + 5;
124 :     int b = floor(pos[1] * 50) + 5;
125 : nseltzer 1608 for (int i = 0; i < 4; i++) {
126 :     for (int j = 0; j < 4; j++) {
127 : nseltzer 1592 tempx[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 0);
128 :     tempy[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 1);
129 :     }
130 :     }
131 :     double t = pos[0] - floor(pos[0]);
132 :     double tempxx[4];
133 :     double tempyx[4];
134 : nseltzer 1608 for (int i = 0; i < 4; i++) {
135 : nseltzer 1592 tempxx[i] = 0.5 * (2 * tempx[0][i] +
136 :     (-tempx[0][i] + tempx[2][i]) * t +
137 :     (2 * tempx[0][i] - 5 * tempx[1][i] + 4 * tempx[2][i] - tempx[3][i]) * t * t +
138 :     (-tempx[0][i] + 3 * tempx[1][i] - 3 * tempx[2][i] + tempx[3][i]) * t * t * t);
139 :     tempyx[i] = 0.5 * (2 * tempy[0][i] * 0 +
140 :     (-tempy[0][i] + tempy[2][i]) * 1 +
141 :     (2 * tempy[0][i] - 5 * tempy[1][i] + 4 * tempy[2][i] - tempy[3][i]) * t * 2 +
142 :     (-tempy[0][i] + 3 * tempy[1][i] - 3 * tempy[2][i] + tempy[3][i]) * t * t * 3);
143 :     }
144 :     t = pos[1] - floor(pos[1]);
145 :     temp[0] = 0.5 * (2 * tempxx[0] * 0 +
146 :     (-tempxx[0] + tempxx[2]) * 1 +
147 :     (2 * tempxx[0] - 5 * tempxx[1] + 4 * tempxx[2] - tempxx[3]) * t * 2 +
148 :     (-tempxx[0] + 3 * tempxx[1] - 3 * tempxx[2] + tempxx[3]) * t * t * 3);
149 :     temp[1] = 0.5 * (2 * tempyx[0] +
150 :     (-tempyx[0] + tempyx[2]) * t +
151 :     (2 * tempyx[0] - 5 * tempyx[1] + 4 * tempyx[2] - tempyx[3]) * t * t +
152 :     (-tempyx[0] + 3 * tempyx[1] - 3 * tempyx[2] + tempyx[3]) * t * t * t);
153 :     return temp[1] - temp[0];
154 :     }
155 :    
156 :     STATIC_INLINE double lerp(double out_min, double out_max, double in_min, double in, double in_max)
157 :     {
158 :     return (out_min * (in_max - in) + out_max * (in - in_min)) / (in_max - in_min);
159 :     }
160 :    
161 :     STATIC_INLINE unsigned char ulerp(double in_min, double in, double in_max)
162 :     {
163 :     return (unsigned char)((0.0f * (in_max - in) + 255.0f * (in - in_min)) / (in_max - in_min));
164 :     }
165 :    
166 : glk 1612 int main (int argc, const char *argv[])
167 : nseltzer 1592 {
168 : glk 1612 const char *me;
169 : nseltzer 1592 char *inFile = "../../data/v.nrrd";
170 :     char *randFile = "../../data/R.nrrd";
171 : glk 1612 char *outFile = "bmark-teem.nrrd";
172 : nseltzer 1592 int imgSizeX = 1020;
173 :     int imgSizeY = 561;
174 :     double h = 0.005;
175 :     int stepNum = 10;
176 :    
177 : glk 1612 me = argv[0];
178 :     airArray *mop;
179 :     mop = airMopNew();
180 :    
181 : nseltzer 1592 // Read in the data
182 :     char* err;
183 : glk 1612 Nrrd *V = nrrdNew();
184 :     airMopAdd(mop, V, (airMopper)nrrdNuke, airMopAlways);
185 :     Nrrd *R = nrrdNew();
186 :     airMopAdd(mop, R, (airMopper)nrrdNuke, airMopAlways);
187 :     if (nrrdLoad(V, inFile, NULL)
188 :     || nrrdLoad(R, randFile, NULL)) {
189 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
190 :     fprintf(stderr, "%s: Trouble reading:\n%s", me, err);
191 :     airMopError(mop);
192 : nseltzer 1592 return -1;
193 :     }
194 :    
195 : glk 1612 Nrrd *nout = nrrdNew();
196 :     airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);
197 :     if (nrrdAlloc_va(nout, nrrdTypeDouble, 3,
198 :     AIR_CAST(size_t, 3),
199 :     AIR_CAST(size_t, 1020),
200 :     AIR_CAST(size_t, 561))) {
201 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
202 :     fprintf(stderr, "%s: Trouble allocating:\n%s", me, err);
203 :     airMopError(mop);
204 : nseltzer 1592 return -1;
205 :     }
206 : glk 1612 double *out_data = AIR_CAST(double *, nout->data);
207 : nseltzer 1592
208 : glk 1601 double t0 = airTime(); // start timing
209 : nseltzer 1592 double pos0[2];
210 :     double forw[2];
211 :     double back[2];
212 :     double temp[2];
213 :     double sum;
214 :     double crl;
215 :     double s, l, c, m;
216 :     int step;
217 :     int num;
218 :     for (int yi = 0; yi < 561; yi++) {
219 :     pos0[1] = lerp(0.04, 3.7, 0, yi, 561);
220 :     for (int xi = 0; xi < 1020; xi++) {
221 :     pos0[0] = lerp(0.04, 6.76, 0, xi, 1020);
222 :     forw[0] = pos0[0];
223 :     forw[1] = pos0[1];
224 :     back[0] = pos0[0];
225 :     back[1] = pos0[1];
226 :     atR(R->data, pos0, &sum);
227 :     sum = -sum;
228 :     step = -1;
229 :     num = -1;
230 :     while (true) {
231 :     if (inside(forw)) {
232 :     atR(R->data, forw, temp);
233 :     sum += temp[0];
234 :     num++;
235 :    
236 :     at(V->data, forw, temp);
237 :     norm2(temp, temp);
238 :     scale2(0.5 * h, temp, temp);
239 :     add2(forw, temp, temp);
240 :     if (inside(temp)) {
241 :     at(V->data, temp, temp);
242 :     norm2(temp, temp);
243 :     scale2(h, temp, temp);
244 :     add2(forw, temp, forw);
245 :     }
246 :     }
247 :     if (inside(back)) {
248 :     atR(R->data, back, temp);
249 :     sum += temp[0];
250 :     num++;
251 :    
252 :     at(V->data, back, temp);
253 :     norm2(temp, temp);
254 :     scale2(-0.5 * h, temp, temp);
255 :     add2(back, temp, temp);
256 :     if (inside(temp)) {
257 :     at(V->data, temp, temp);
258 :     norm2(temp, temp);
259 :     scale2(-1 * h, temp, temp);
260 :     add2(back, temp, back);
261 :     }
262 :     }
263 :     step += 1;
264 :     if ((!inside(forw) && !inside(back)) || step == stepNum) {
265 :     at(V->data, pos0, temp);
266 :     sum = pow(mag2(temp), 0.3) * sum / num;
267 :     crl = crl_at(V->data, pos0);
268 :     s = lerp(0, 1, 0, (crl>0)?crl:-1*crl, 2.794);
269 : nseltzer 1611 l = lerp(0.0, 1.0, -1.6, sum, 2.6);
270 : nseltzer 1592
271 :     c = (1 - (((2 * l - 1)>0)?(2 * l - 1):(1 - 2 * l))) * s;
272 :     m = l - c / 2;
273 : glk 1612 if (crl > 0) {
274 :     ELL_3V_SET(out_data + 3*(xi + 1020*yi), c+m, m, c+m);
275 : nseltzer 1592 }
276 :     else {
277 : glk 1612 ELL_3V_SET(out_data + 3*(xi + 1020*yi), m, c+m, m);
278 : nseltzer 1592 }
279 :     break;
280 :     }
281 :     }
282 :     }
283 :     }
284 :    
285 : glk 1601 double totalTime = airTime() - t0;
286 : nseltzer 1592 printf("usr=%f\n", totalTime);
287 :    
288 : glk 1612 if (nrrdSave(outFile, nout, NULL)) {
289 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
290 :     fprintf(stderr, "%s: Trouble saving:\n%s", me, err);
291 :     airMopError(mop);
292 : nseltzer 1592 return -1;
293 :     }
294 :    
295 : glk 1612 airMopOkay(mop);
296 : nseltzer 1592 return 0;
297 :     }

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