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

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