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 3349 - (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 : jhr 3349 * This code is part of the Diderot Project (http://diderot-language.cs.uchicago.edu)
8 :     *
9 :     * COPYRIGHT (c) 2015 The University of Chicago
10 : nseltzer 1592 * All rights reserved.
11 :     */
12 :    
13 :     #include <stdio.h>
14 :     #include <math.h>
15 :     #include <stdbool.h>
16 :     #include "teem/nrrd.h"
17 : glk 1612 #include "teem/ell.h"
18 : nseltzer 1592 #include "teem/gage.h"
19 :    
20 : jhr 1596 #include "teem-defs.h"
21 : nseltzer 1592
22 :     STATIC_INLINE void scale2(double scl, double vec[2], double res[2])
23 :     {
24 :     res[0] = scl * vec[0];
25 :     res[1] = scl * vec[1];
26 :     }
27 :    
28 :     STATIC_INLINE void add2(double vec1[2], double vec2[2], double res[2])
29 :     {
30 :     res[0] = vec1[0] + vec2[0];
31 :     res[1] = vec1[1] + vec2[1];
32 :     }
33 :    
34 :     STATIC_INLINE double dot2(double vec1[2], double vec2[2])
35 :     {
36 :     return (vec1[0] * vec2[0] + vec1[1] * vec2[1]);
37 :     }
38 :    
39 :     STATIC_INLINE double mag2(double vec[2])
40 :     {
41 :     return sqrt(dot2(vec, vec));
42 :     }
43 :    
44 :     STATIC_INLINE void norm2(double vec[2], double res[2])
45 :     {
46 :     double mag = mag2(vec);
47 :     res[0] = vec[0] / mag;
48 :     res[1] = vec[1] / mag;
49 :     }
50 :    
51 :     STATIC_INLINE bool inside(double pos[2])
52 :     {
53 :     // XXX - Hack
54 : nseltzer 1614 if (pos[0] < 0.00 || pos[0] >= 6.78) return false;
55 :     if (pos[1] < 0.00 || pos[1] >= 3.72) return false;
56 : nseltzer 1592 return true;
57 :     }
58 :    
59 :     void atR(void *_data, double pos[2], double *out)
60 :     {
61 :     float *data = (float*)_data;
62 : nseltzer 1608 double temp[2][2];
63 : nseltzer 1614 int a = floor(pos[0] * 150) + 1;
64 :     int b = floor(pos[1] * 150) + 1;
65 :     for (int i = 0; i < 2; i++) {
66 : nseltzer 1608 for (int j = 0; j < 2; j++) {
67 :     temp[i][j] = *(data + (i + a) + (j + b) * 1020);
68 :     }
69 :     }
70 :     double temp2[2];
71 :     double t = pos[0]*150 - floor(pos[0]*150);
72 :     for (int i = 0; i < 2; i++) {
73 :     temp2[i] = (temp[0][i] * (1 - t) + temp[1][i] * t);
74 :     }
75 :     t = pos[1]*150 - floor(pos[1]*150);
76 :     *out = (temp2[0] * (1 - t) + temp2[1] * t);
77 : nseltzer 1592 }
78 :    
79 :     void at(void* _data, double pos[2], double out[2])
80 :     {
81 :     // XXX - Hack
82 :     float *data = (float*)_data;
83 :     double tempx[4][4];
84 :     double tempy[4][4];
85 :     int a = floor(pos[0] * 50) + 5;
86 :     int b = floor(pos[1] * 50) + 5;
87 : nseltzer 1608 for (int i = 0; i < 4; i++) {
88 :     for (int j = 0; j < 4; j++) {
89 : nseltzer 1592 tempx[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 0);
90 :     tempy[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 1);
91 :     }
92 :     }
93 : nseltzer 1614 double t = pos[0]*50 - floor(pos[0]*50);
94 : nseltzer 1592 double tempxx[4];
95 :     double tempyx[4];
96 : nseltzer 1608 for (int i = 0; i < 4; i++) {
97 : nseltzer 1614 tempxx[i] = (tempx[0][i] + 4 * tempx[1][i] + tempx[2][i] +
98 :     (-3 * tempx[0][i] + 3 * tempx[2][i]) * t +
99 :     (3 * tempx[0][i] - 6 * tempx[1][i] + 3 * tempx[2][i]) * t * t +
100 :     (-1 * tempx[0][i] + 3 * tempx[1][i] - 3 * tempx[2][i] + tempx[3][i]) * t * t * t) / 6.0;
101 :     tempyx[i] = (tempy[0][i] + 4 * tempy[1][i] + tempy[2][i] +
102 :     (-3 * tempy[0][i] + 3 * tempy[2][i]) * t +
103 :     (3 * tempy[0][i] - 6 * tempy[1][i] + 3 * tempy[2][i]) * t * t +
104 :     (-1 * tempy[0][i] + 3 * tempy[1][i] - 3 * tempy[2][i] + tempy[3][i]) * t * t * t) / 6.0;
105 : nseltzer 1592 }
106 : nseltzer 1614 t = pos[1]*50 - floor(pos[1]*50);
107 :     out[0] = (tempxx[0] + 4 * tempxx[1] + tempxx[2] +
108 :     (-3 * tempxx[0] + 3 * tempxx[2]) * t +
109 :     (3 * tempxx[0] - 6 * tempxx[1] + 3 * tempxx[2]) * t * t +
110 :     (-1 * tempxx[0] + 3 * tempxx[1] - 3 * tempxx[2] + tempxx[3]) * t * t * t) / 6.0;
111 :     out[1] = (tempyx[0] + 4 * tempyx[1] + tempyx[2] +
112 :     (-3 * tempyx[0] + 3 * tempyx[2]) * t +
113 :     (3 * tempyx[0] - 6 * tempyx[1] + 3 * tempyx[2]) * t * t +
114 :     (-1 * tempyx[0] + 3 * tempyx[1] - 3 * tempyx[2] + tempyx[3]) * t * t * t) / 6.0;
115 : nseltzer 1592 return;
116 :     }
117 :    
118 :     double crl_at(void* _data, double pos[2])
119 :     {
120 :     // XXX - Hack
121 :     float *data = (float*)_data;
122 :     double tempx[4][4];
123 :     double tempy[4][4];
124 :     double temp[2];
125 :     int a = floor(pos[0] * 50) + 5;
126 :     int b = floor(pos[1] * 50) + 5;
127 : nseltzer 1608 for (int i = 0; i < 4; i++) {
128 :     for (int j = 0; j < 4; j++) {
129 : nseltzer 1592 tempx[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 0);
130 :     tempy[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 1);
131 :     }
132 :     }
133 : nseltzer 1614 double t = pos[0]*50 - floor(pos[0]*50);
134 : nseltzer 1592 double tempxx[4];
135 :     double tempyx[4];
136 : nseltzer 1608 for (int i = 0; i < 4; i++) {
137 : nseltzer 1614 tempxx[i] = (tempx[0][i] + 4 * tempx[1][i] + tempx[2][i] +
138 :     (-3 * tempx[0][i] + 3 * tempx[2][i]) * t +
139 :     (3 * tempx[0][i] - 6 * tempx[1][i] + 3 * tempx[2][i]) * t * t +
140 :     (-1 * tempx[0][i] + 3 * tempx[1][i] - 3 * tempx[2][i] + tempx[3][i]) * t * t * t) / 6.0;
141 :     tempyx[i] = ((-3 * tempy[0][i] + 3 * tempy[2][i]) * 1 +
142 :     (3 * tempy[0][i] - 6 * tempy[1][i] + 3 * tempy[2][i]) * t * 2 +
143 :     (-1 * tempy[0][i] + 3 * tempy[1][i] - 3 * tempy[2][i] + tempy[3][i]) * t * t * 3) / 6.0;
144 : nseltzer 1592 }
145 : nseltzer 1614 t = pos[1]*50 - floor(pos[1]*50);
146 :     temp[0] = ((-3 * tempxx[0] + 3 * tempxx[2]) * 1 +
147 :     (3 * tempxx[0] - 6 * tempxx[1] + 3 * tempxx[2]) * t * 2 +
148 :     (-1 * tempxx[0] + 3 * tempxx[1] - 3 * tempxx[2] + tempxx[3]) * t * t * 3) / 6.0;
149 :     temp[1] = (tempyx[0] + 4 * tempyx[1] + tempyx[2] +
150 :     (-3 * tempyx[0] + 3 * tempyx[2]) * t +
151 :     (3 * tempyx[0] - 6 * tempyx[1] + 3 * tempyx[2]) * t * t +
152 :     (-1 * tempyx[0] + 3 * tempyx[1] - 3 * tempyx[2] + tempyx[3]) * t * t * t) / 6.0;
153 : nseltzer 1592 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 : nseltzer 1614 pos0[1] = lerp(0.04, 3.7, 0, yi, 560);
220 : nseltzer 1592 for (int xi = 0; xi < 1020; xi++) {
221 : nseltzer 1614 pos0[0] = lerp(0.04, 6.76, 0, xi, 1019);
222 : nseltzer 1592 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 : nseltzer 1614 s = lerp(0, 1, 0, (crl>0)?crl:-1*crl, 2.64);
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