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

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