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 1608 - (view) (download) (as text)

1 : nseltzer 1592 /*! \file /lic2d/bmark-teem.c
2 :     *
3 :     * \author Nick Seltzer
4 :     */
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 :     #include "teem/gage.h"
16 :    
17 : jhr 1596 #include "teem-defs.h"
18 : nseltzer 1592
19 :     STATIC_INLINE void scale2(double scl, double vec[2], double res[2])
20 :     {
21 :     res[0] = scl * vec[0];
22 :     res[1] = scl * vec[1];
23 :     }
24 :    
25 :     STATIC_INLINE void add2(double vec1[2], double vec2[2], double res[2])
26 :     {
27 :     res[0] = vec1[0] + vec2[0];
28 :     res[1] = vec1[1] + vec2[1];
29 :     }
30 :    
31 :     STATIC_INLINE double dot2(double vec1[2], double vec2[2])
32 :     {
33 :     return (vec1[0] * vec2[0] + vec1[1] * vec2[1]);
34 :     }
35 :    
36 :     STATIC_INLINE double mag2(double vec[2])
37 :     {
38 :     return sqrt(dot2(vec, vec));
39 :     }
40 :    
41 :     STATIC_INLINE void norm2(double vec[2], double res[2])
42 :     {
43 :     double mag = mag2(vec);
44 :     res[0] = vec[0] / mag;
45 :     res[1] = vec[1] / mag;
46 :     }
47 :    
48 :     STATIC_INLINE bool inside(double pos[2])
49 :     {
50 :     // XXX - Hack
51 : nseltzer 1608 if (pos[0] < 0.04 || pos[0] > 6.76) return false;
52 :     if (pos[1] < 0.04 || pos[1] > 3.7) return false;
53 : nseltzer 1592 return true;
54 :     }
55 :    
56 :     void atR(void *_data, double pos[2], double *out)
57 :     {
58 :     float *data = (float*)_data;
59 : nseltzer 1608 double temp[2][2];
60 : nseltzer 1592 int a = floor(pos[0] * 150);
61 :     int b = floor(pos[1] * 150);
62 : nseltzer 1608 for (int i = 0; i < 4; i++) {
63 :     for (int j = 0; j < 2; j++) {
64 :     temp[i][j] = *(data + (i + a) + (j + b) * 1020);
65 :     }
66 :     }
67 :     double temp2[2];
68 :     double t = pos[0]*150 - floor(pos[0]*150);
69 :     for (int i = 0; i < 2; i++) {
70 :     temp2[i] = (temp[0][i] * (1 - t) + temp[1][i] * t);
71 :     }
72 :     t = pos[1]*150 - floor(pos[1]*150);
73 :     *out = (temp2[0] * (1 - t) + temp2[1] * t);
74 : nseltzer 1592 }
75 :    
76 :     void at(void* _data, double pos[2], double out[2])
77 :     {
78 :     // XXX - Hack
79 :     float *data = (float*)_data;
80 :     double tempx[4][4];
81 :     double tempy[4][4];
82 :     int a = floor(pos[0] * 50) + 5;
83 :     int b = floor(pos[1] * 50) + 5;
84 : nseltzer 1608 for (int i = 0; i < 4; i++) {
85 :     for (int j = 0; j < 4; j++) {
86 : nseltzer 1592 tempx[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 0);
87 :     tempy[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 1);
88 :     }
89 :     }
90 : nseltzer 1608 double t = pos[0]*150 - floor(pos[0]*150);
91 : nseltzer 1592 double tempxx[4];
92 :     double tempyx[4];
93 : nseltzer 1608 for (int i = 0; i < 4; i++) {
94 : nseltzer 1592 tempxx[i] = 0.5 * (2 * tempx[0][i] +
95 :     (-tempx[0][i] + tempx[2][i]) * t +
96 :     (2 * tempx[0][i] - 5 * tempx[1][i] + 4 * tempx[2][i] - tempx[3][i]) * t * t +
97 :     (-tempx[0][i] + 3 * tempx[1][i] - 3 * tempx[2][i] + tempx[3][i]) * t * t * t);
98 :     tempyx[i] = 0.5 * (2 * tempy[0][i] +
99 :     (-tempy[0][i] + tempy[2][i]) * t +
100 :     (2 * tempy[0][i] - 5 * tempy[1][i] + 4 * tempy[2][i] - tempy[3][i]) * t * t +
101 :     (-tempy[0][i] + 3 * tempy[1][i] - 3 * tempy[2][i] + tempy[3][i]) * t * t * t);
102 :     }
103 : nseltzer 1608 t = pos[1]*150 - floor(pos[1]*150);
104 : nseltzer 1592 out[0] = 0.5 * (2 * tempxx[0] +
105 :     (-tempxx[0] + tempxx[2]) * t +
106 :     (2 * tempxx[0] - 5 * tempxx[1] + 4 * tempxx[2] - tempxx[3]) * t * t +
107 :     (-tempxx[0] + 3 * tempxx[1] - 3 * tempxx[2] + tempxx[3]) * t * t * t);
108 :     out[1] = 0.5 * (2 * tempyx[0] +
109 :     (-tempyx[0] + tempyx[2]) * t +
110 :     (2 * tempyx[0] - 5 * tempyx[1] + 4 * tempyx[2] - tempyx[3]) * t * t +
111 :     (-tempyx[0] + 3 * tempyx[1] - 3 * tempyx[2] + tempyx[3]) * t * t * t);
112 :     return;
113 :     }
114 :    
115 :     double crl_at(void* _data, double pos[2])
116 :     {
117 :     // XXX - Hack
118 :     float *data = (float*)_data;
119 :     double tempx[4][4];
120 :     double tempy[4][4];
121 :     double temp[2];
122 :     int a = floor(pos[0] * 50) + 5;
123 :     int b = floor(pos[1] * 50) + 5;
124 : nseltzer 1608 for (int i = 0; i < 4; i++) {
125 :     for (int j = 0; j < 4; j++) {
126 : nseltzer 1592 tempx[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 0);
127 :     tempy[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 1);
128 :     }
129 :     }
130 :     double t = pos[0] - floor(pos[0]);
131 :     double tempxx[4];
132 :     double tempyx[4];
133 : nseltzer 1608 for (int i = 0; i < 4; i++) {
134 : nseltzer 1592 tempxx[i] = 0.5 * (2 * tempx[0][i] +
135 :     (-tempx[0][i] + tempx[2][i]) * t +
136 :     (2 * tempx[0][i] - 5 * tempx[1][i] + 4 * tempx[2][i] - tempx[3][i]) * t * t +
137 :     (-tempx[0][i] + 3 * tempx[1][i] - 3 * tempx[2][i] + tempx[3][i]) * t * t * t);
138 :     tempyx[i] = 0.5 * (2 * tempy[0][i] * 0 +
139 :     (-tempy[0][i] + tempy[2][i]) * 1 +
140 :     (2 * tempy[0][i] - 5 * tempy[1][i] + 4 * tempy[2][i] - tempy[3][i]) * t * 2 +
141 :     (-tempy[0][i] + 3 * tempy[1][i] - 3 * tempy[2][i] + tempy[3][i]) * t * t * 3);
142 :     }
143 :     t = pos[1] - floor(pos[1]);
144 :     temp[0] = 0.5 * (2 * tempxx[0] * 0 +
145 :     (-tempxx[0] + tempxx[2]) * 1 +
146 :     (2 * tempxx[0] - 5 * tempxx[1] + 4 * tempxx[2] - tempxx[3]) * t * 2 +
147 :     (-tempxx[0] + 3 * tempxx[1] - 3 * tempxx[2] + tempxx[3]) * t * t * 3);
148 :     temp[1] = 0.5 * (2 * tempyx[0] +
149 :     (-tempyx[0] + tempyx[2]) * t +
150 :     (2 * tempyx[0] - 5 * tempyx[1] + 4 * tempyx[2] - tempyx[3]) * t * t +
151 :     (-tempyx[0] + 3 * tempyx[1] - 3 * tempyx[2] + tempyx[3]) * t * t * t);
152 :     return temp[1] - temp[0];
153 :     }
154 :    
155 :     STATIC_INLINE double lerp(double out_min, double out_max, double in_min, double in, double in_max)
156 :     {
157 :     return (out_min * (in_max - in) + out_max * (in - in_min)) / (in_max - in_min);
158 :     }
159 :    
160 :     STATIC_INLINE unsigned char ulerp(double in_min, double in, double in_max)
161 :     {
162 :     return (unsigned char)((0.0f * (in_max - in) + 255.0f * (in - in_min)) / (in_max - in_min));
163 :     }
164 :    
165 :     int main ()
166 :     {
167 :     char *inFile = "../../data/v.nrrd";
168 :     char *randFile = "../../data/R.nrrd";
169 :     char *outFile = "lic2d.ppm";
170 :     int imgSizeX = 1020;
171 :     int imgSizeY = 561;
172 :     double h = 0.005;
173 :     int stepNum = 10;
174 :    
175 :     // Read in the data
176 :     char* err;
177 :     Nrrd *V;
178 :     int status;
179 :     V = nrrdNew();
180 : nseltzer 1608 if (V == NULL) {
181 : nseltzer 1592 err = biffGetDone(NRRD);
182 :     fprintf(stderr, "Trouble allocating nrrd struct:\n%s", err);
183 :     free(err);
184 :     return -1;
185 :     }
186 :    
187 :     status = nrrdLoad(V, inFile, NULL);
188 :     if (status) {
189 :     err = biffGetDone(NRRD);
190 :     fprintf(stderr, "Trouble reading \"%s\":\n%s", inFile, err);
191 :     free(err);
192 :     nrrdNix(V);
193 :     return -1;
194 :     }
195 :    
196 :     // Read in the noise
197 :     Nrrd *R;
198 :     R = nrrdNew();
199 : nseltzer 1608 if (R == NULL) {
200 : nseltzer 1592 err = biffGetDone(NRRD);
201 :     fprintf(stderr, "Trouble allocating nrrd struct:\n%s", err);
202 :     free(err);
203 :     nrrdNuke(V);
204 :     return -1;
205 :     }
206 :    
207 :     status = nrrdLoad(R, randFile, NULL);
208 :     if (status) {
209 :     err = biffGetDone(NRRD);
210 :     fprintf(stderr, "Trouble reading \"%s\":\n%s", randFile, err);
211 :     free(err);
212 :     nrrdNix(R);
213 :     nrrdNuke(V);
214 :     return -1;
215 :     }
216 :    
217 :     double *out_data = malloc(sizeof(double) * 3 * 1020 * 561);
218 :     if (out_data == NULL) {
219 :     fprintf(stderr, "Trouble with malloc\n");
220 :     nrrdNuke(V);
221 :     nrrdNuke(R);
222 :     return -1;
223 :     }
224 :    
225 : glk 1601 double t0 = airTime(); // start timing
226 : nseltzer 1592 double pos0[2];
227 :     double forw[2];
228 :     double back[2];
229 :     double temp[2];
230 :     double sum;
231 :     double crl;
232 :     double s, l, c, m;
233 :     int step;
234 :     int num;
235 :     for (int yi = 0; yi < 561; yi++) {
236 :     pos0[1] = lerp(0.04, 3.7, 0, yi, 561);
237 :     for (int xi = 0; xi < 1020; xi++) {
238 :     pos0[0] = lerp(0.04, 6.76, 0, xi, 1020);
239 :     forw[0] = pos0[0];
240 :     forw[1] = pos0[1];
241 :     back[0] = pos0[0];
242 :     back[1] = pos0[1];
243 :     atR(R->data, pos0, &sum);
244 :     sum = -sum;
245 :     step = -1;
246 :     num = -1;
247 :     while (true) {
248 :     if (inside(forw)) {
249 :     atR(R->data, forw, temp);
250 :     sum += temp[0];
251 :     num++;
252 :    
253 :     at(V->data, forw, temp);
254 :     norm2(temp, temp);
255 :     scale2(0.5 * h, temp, temp);
256 :     add2(forw, temp, temp);
257 :     if (inside(temp)) {
258 :     at(V->data, temp, temp);
259 :     norm2(temp, temp);
260 :     scale2(h, temp, temp);
261 :     add2(forw, temp, forw);
262 :     }
263 :     }
264 :     if (inside(back)) {
265 :     atR(R->data, back, temp);
266 :     sum += temp[0];
267 :     num++;
268 :    
269 :     at(V->data, back, temp);
270 :     norm2(temp, temp);
271 :     scale2(-0.5 * h, temp, temp);
272 :     add2(back, temp, temp);
273 :     if (inside(temp)) {
274 :     at(V->data, temp, temp);
275 :     norm2(temp, temp);
276 :     scale2(-1 * h, temp, temp);
277 :     add2(back, temp, back);
278 :     }
279 :     }
280 :     step += 1;
281 :     if ((!inside(forw) && !inside(back)) || step == stepNum) {
282 :     at(V->data, pos0, temp);
283 :     sum = pow(mag2(temp), 0.3) * sum / num;
284 :     crl = crl_at(V->data, pos0);
285 :     s = lerp(0, 1, 0, (crl>0)?crl:-1*crl, 2.794);
286 :     l = lerp(0.0, 1.0, -2.583, sum, 2.583);
287 :    
288 :     c = (1 - (((2 * l - 1)>0)?(2 * l - 1):(1 - 2 * l))) * s;
289 :     m = l - c / 2;
290 :     if (crl < 0) {
291 :     *(out_data + xi * 3 + yi * 3060 + 0) = c + m;
292 :     *(out_data + xi * 3 + yi * 3060 + 1) = m;
293 :     *(out_data + xi * 3 + yi * 3060 + 2) = c + m;
294 :     }
295 :     else {
296 :     *(out_data + xi * 3 + yi * 3060 + 0) = m;
297 :     *(out_data + xi * 3 + yi * 3060 + 1) = c + m;
298 :     *(out_data + xi * 3 + yi * 3060 + 2) = m;
299 :     }
300 :     break;
301 :     }
302 :     }
303 :     }
304 :     }
305 :    
306 :     nrrdNuke(V);
307 :     nrrdNuke(R);
308 :    
309 : glk 1601 double totalTime = airTime() - t0;
310 : nseltzer 1592 printf("usr=%f\n", totalTime);
311 :    
312 :     unsigned char *uc_out_data = malloc(sizeof(unsigned char) * 3 * 1020 * 561);
313 :     if (uc_out_data == NULL) {
314 :     fprintf(stderr, "Trouble with malloc\n");
315 :     free(out_data);
316 :     }
317 :     for (int i = 0; i < 3; i++) {
318 :     for (int j = 0; j < 1020; j++) {
319 :     for (int k = 0; k < 561; k++) {
320 :     *(uc_out_data + k * 3060 + j * 3 + i) = ulerp(0, *(out_data + k * 3060 + j * 3 + i), 1);
321 :     }
322 :     }
323 :     }
324 :    
325 :     free(out_data);
326 :    
327 :     // Write out data
328 :     Nrrd *nout;
329 :     nout = nrrdNew();
330 :     if (nout == NULL) {
331 :     err = biffGetDone(NRRD);
332 :     fprintf(stderr, "Trouble allocating nrrd struct:\n%s", err);
333 :     free(err);
334 :     free(uc_out_data);
335 :     return -1;
336 :     }
337 :    
338 :     size_t sizes[3] = {3, 1020, 561};
339 :     status = nrrdWrap_nva(nout, uc_out_data, nrrdTypeUChar, 3, sizes);
340 :     if (status) {
341 :     err = biffGetDone(NRRD);
342 :     fprintf(stderr, "Trouble wrapping nrrd struct:\n%s", err);
343 :     free(err);
344 :     free(uc_out_data);
345 :     nrrdNix(nout);
346 :     return -1;
347 :     }
348 :    
349 :     status = nrrdSave(outFile, nout, NULL);
350 :     if (status) {
351 :     err = biffGetDone(NRRD);
352 :     fprintf(stderr, "Trouble wrapping nrrd struct:\n%s", err);
353 :     free(err);
354 :     nrrdNuke(nout);
355 :     return -1;
356 :     }
357 :    
358 :     nrrdNuke(nout);
359 :     return 0;
360 :     }

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