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 1592 - (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 :     #include <sys/time.h>
18 :    
19 :     static double GetTime ()
20 :     {
21 :     struct timeval tv;
22 :    
23 :     gettimeofday (&tv, 0);
24 :    
25 :     return (double)tv.tv_sec + 0.000001 * (double)tv.tv_usec;
26 :     }
27 :    
28 :     #define STATIC_INLINE static inline
29 :    
30 :     STATIC_INLINE void scale2(double scl, double vec[2], double res[2])
31 :     {
32 :     res[0] = scl * vec[0];
33 :     res[1] = scl * vec[1];
34 :     }
35 :    
36 :     STATIC_INLINE void add2(double vec1[2], double vec2[2], double res[2])
37 :     {
38 :     res[0] = vec1[0] + vec2[0];
39 :     res[1] = vec1[1] + vec2[1];
40 :     }
41 :    
42 :     STATIC_INLINE double dot2(double vec1[2], double vec2[2])
43 :     {
44 :     return (vec1[0] * vec2[0] + vec1[1] * vec2[1]);
45 :     }
46 :    
47 :     STATIC_INLINE double mag2(double vec[2])
48 :     {
49 :     return sqrt(dot2(vec, vec));
50 :     }
51 :    
52 :     STATIC_INLINE void norm2(double vec[2], double res[2])
53 :     {
54 :     double mag = mag2(vec);
55 :     res[0] = vec[0] / mag;
56 :     res[1] = vec[1] / mag;
57 :     }
58 :    
59 :     STATIC_INLINE bool inside(double pos[2])
60 :     {
61 :     // XXX - Hack
62 :     if(pos[0] < 0.04 || pos[0] > 6.76) return false;
63 :     if(pos[1] < 0.04 || pos[1] > 3.7) return false;
64 :     return true;
65 :     }
66 :    
67 :     void atR(void *_data, double pos[2], double *out)
68 :     {
69 :     float *data = (float*)_data;
70 :     int a = floor(pos[0] * 150);
71 :     int b = floor(pos[1] * 150);
72 :     *out = *(data + a + b * 1020);
73 :     }
74 :    
75 :     void at(void* _data, double pos[2], double out[2])
76 :     {
77 :     // XXX - Hack
78 :     float *data = (float*)_data;
79 :     double tempx[4][4];
80 :     double tempy[4][4];
81 :     int a = floor(pos[0] * 50) + 5;
82 :     int b = floor(pos[1] * 50) + 5;
83 :     for(int i = 0; i < 4; i++) {
84 :     for(int j = 0; j < 4; j++) {
85 :     tempx[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 0);
86 :     tempy[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 1);
87 :     }
88 :     }
89 :     double t = pos[0] - floor(pos[0]);
90 :     double tempxx[4];
91 :     double tempyx[4];
92 :     for(int i = 0; i < 4; i++) {
93 :     tempxx[i] = 0.5 * (2 * tempx[0][i] +
94 :     (-tempx[0][i] + tempx[2][i]) * t +
95 :     (2 * tempx[0][i] - 5 * tempx[1][i] + 4 * tempx[2][i] - tempx[3][i]) * t * t +
96 :     (-tempx[0][i] + 3 * tempx[1][i] - 3 * tempx[2][i] + tempx[3][i]) * t * t * t);
97 :     tempyx[i] = 0.5 * (2 * tempy[0][i] +
98 :     (-tempy[0][i] + tempy[2][i]) * t +
99 :     (2 * tempy[0][i] - 5 * tempy[1][i] + 4 * tempy[2][i] - tempy[3][i]) * t * t +
100 :     (-tempy[0][i] + 3 * tempy[1][i] - 3 * tempy[2][i] + tempy[3][i]) * t * t * t);
101 :     }
102 :     t = pos[1] - floor(pos[1]);
103 :     out[0] = 0.5 * (2 * tempxx[0] +
104 :     (-tempxx[0] + tempxx[2]) * t +
105 :     (2 * tempxx[0] - 5 * tempxx[1] + 4 * tempxx[2] - tempxx[3]) * t * t +
106 :     (-tempxx[0] + 3 * tempxx[1] - 3 * tempxx[2] + tempxx[3]) * t * t * t);
107 :     out[1] = 0.5 * (2 * tempyx[0] +
108 :     (-tempyx[0] + tempyx[2]) * t +
109 :     (2 * tempyx[0] - 5 * tempyx[1] + 4 * tempyx[2] - tempyx[3]) * t * t +
110 :     (-tempyx[0] + 3 * tempyx[1] - 3 * tempyx[2] + tempyx[3]) * t * t * t);
111 :     return;
112 :     }
113 :    
114 :     double crl_at(void* _data, double pos[2])
115 :     {
116 :     // XXX - Hack
117 :     float *data = (float*)_data;
118 :     double tempx[4][4];
119 :     double tempy[4][4];
120 :     double temp[2];
121 :     int a = floor(pos[0] * 50) + 5;
122 :     int b = floor(pos[1] * 50) + 5;
123 :     for(int i = 0; i < 4; i++) {
124 :     for(int j = 0; j < 4; j++) {
125 :     tempx[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 0);
126 :     tempy[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 1);
127 :     }
128 :     }
129 :     double t = pos[0] - floor(pos[0]);
130 :     double tempxx[4];
131 :     double tempyx[4];
132 :     for(int i = 0; i < 4; i++) {
133 :     tempxx[i] = 0.5 * (2 * tempx[0][i] +
134 :     (-tempx[0][i] + tempx[2][i]) * t +
135 :     (2 * tempx[0][i] - 5 * tempx[1][i] + 4 * tempx[2][i] - tempx[3][i]) * t * t +
136 :     (-tempx[0][i] + 3 * tempx[1][i] - 3 * tempx[2][i] + tempx[3][i]) * t * t * t);
137 :     tempyx[i] = 0.5 * (2 * tempy[0][i] * 0 +
138 :     (-tempy[0][i] + tempy[2][i]) * 1 +
139 :     (2 * tempy[0][i] - 5 * tempy[1][i] + 4 * tempy[2][i] - tempy[3][i]) * t * 2 +
140 :     (-tempy[0][i] + 3 * tempy[1][i] - 3 * tempy[2][i] + tempy[3][i]) * t * t * 3);
141 :     }
142 :     t = pos[1] - floor(pos[1]);
143 :     temp[0] = 0.5 * (2 * tempxx[0] * 0 +
144 :     (-tempxx[0] + tempxx[2]) * 1 +
145 :     (2 * tempxx[0] - 5 * tempxx[1] + 4 * tempxx[2] - tempxx[3]) * t * 2 +
146 :     (-tempxx[0] + 3 * tempxx[1] - 3 * tempxx[2] + tempxx[3]) * t * t * 3);
147 :     temp[1] = 0.5 * (2 * tempyx[0] +
148 :     (-tempyx[0] + tempyx[2]) * t +
149 :     (2 * tempyx[0] - 5 * tempyx[1] + 4 * tempyx[2] - tempyx[3]) * t * t +
150 :     (-tempyx[0] + 3 * tempyx[1] - 3 * tempyx[2] + tempyx[3]) * t * t * t);
151 :     return temp[1] - temp[0];
152 :     }
153 :    
154 :     STATIC_INLINE double lerp(double out_min, double out_max, double in_min, double in, double in_max)
155 :     {
156 :     return (out_min * (in_max - in) + out_max * (in - in_min)) / (in_max - in_min);
157 :     }
158 :    
159 :     STATIC_INLINE unsigned char ulerp(double in_min, double in, double in_max)
160 :     {
161 :     return (unsigned char)((0.0f * (in_max - in) + 255.0f * (in - in_min)) / (in_max - in_min));
162 :     }
163 :    
164 :     int main ()
165 :     {
166 :     char *inFile = "../../data/v.nrrd";
167 :     char *randFile = "../../data/R.nrrd";
168 :     char *outFile = "lic2d.ppm";
169 :     int imgSizeX = 1020;
170 :     int imgSizeY = 561;
171 :     double h = 0.005;
172 :     int stepNum = 10;
173 :    
174 :     // Read in the data
175 :     char* err;
176 :     Nrrd *V;
177 :     int status;
178 :     V = nrrdNew();
179 :     if(V == NULL) {
180 :     err = biffGetDone(NRRD);
181 :     fprintf(stderr, "Trouble allocating nrrd struct:\n%s", err);
182 :     free(err);
183 :     return -1;
184 :     }
185 :    
186 :     status = nrrdLoad(V, inFile, NULL);
187 :     if (status) {
188 :     err = biffGetDone(NRRD);
189 :     fprintf(stderr, "Trouble reading \"%s\":\n%s", inFile, err);
190 :     free(err);
191 :     nrrdNix(V);
192 :     return -1;
193 :     }
194 :    
195 :     // Read in the noise
196 :     Nrrd *R;
197 :     R = nrrdNew();
198 :     if(R == NULL) {
199 :     err = biffGetDone(NRRD);
200 :     fprintf(stderr, "Trouble allocating nrrd struct:\n%s", err);
201 :     free(err);
202 :     nrrdNuke(V);
203 :     return -1;
204 :     }
205 :    
206 :     status = nrrdLoad(R, randFile, NULL);
207 :     if (status) {
208 :     err = biffGetDone(NRRD);
209 :     fprintf(stderr, "Trouble reading \"%s\":\n%s", randFile, err);
210 :     free(err);
211 :     nrrdNix(R);
212 :     nrrdNuke(V);
213 :     return -1;
214 :     }
215 :    
216 :     double *out_data = malloc(sizeof(double) * 3 * 1020 * 561);
217 :     if (out_data == NULL) {
218 :     fprintf(stderr, "Trouble with malloc\n");
219 :     nrrdNuke(V);
220 :     nrrdNuke(R);
221 :     return -1;
222 :     }
223 :    
224 :    
225 :     double t0 = GetTime(); // start timing
226 :     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 :     double totalTime = GetTime() - t0;
310 :     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 :     //status = nrrdWrap_va(nout, uc_out_data, nrrdTypeUChar, 3, 3, 1020, 561);
341 :     // XXX - Teem bug makes nrrdWrap_va not work??
342 :     if (status) {
343 :     err = biffGetDone(NRRD);
344 :     fprintf(stderr, "Trouble wrapping nrrd struct:\n%s", err);
345 :     free(err);
346 :     free(uc_out_data);
347 :     nrrdNix(nout);
348 :     return -1;
349 :     }
350 :    
351 :     status = nrrdSave(outFile, nout, NULL);
352 :     if (status) {
353 :     err = biffGetDone(NRRD);
354 :     fprintf(stderr, "Trouble wrapping nrrd struct:\n%s", err);
355 :     free(err);
356 :     nrrdNuke(nout);
357 :     return -1;
358 :     }
359 :    
360 :     nrrdNuke(nout);
361 :     return 0;
362 :     }

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