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

Parent Directory | 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