Annotation of /branches/pure-cfg/src/lib/diderot.c

Parent Directory | Revision Log

---

Revision 571 - (view) (download) (as text)

1 :	jhr	439	/*! \file diderot.c
2 :			*
3 :			* \author John Reppy
4 :			*/
5 :
6 :			/*
7 :			* COPYRIGHT (c) 2010 The Diderot Project (http://diderot-language.cs.uchicago.edu)
8 :			* All rights reserved.
9 :			*/
10 :
11 :	jhr	566	#include "Diderot/diderot.h"
12 :	jhr	439	#include <teem/nrrd.h>
13 :
14 :			/* Transformation matrix operations */
15 :
16 :			STATIC_INLINE Diderot_real_t DetM3x3 (
17 :			Diderot_real_t a, Diderot_real_t b, Diderot_real_t c,
18 :			Diderot_real_t d, Diderot_real_t e, Diderot_real_t f,
19 :			Diderot_real_t g, Diderot_real_t h, Diderot_real_t i)
20 :			{
21 :			return ( (a)*(e)*(i)
22 :			+ (d)*(h)*(c)
23 :			+ (g)*(b)*(f)
24 :			- (g)*(e)*(c)
25 :			- (d)*(b)*(i)
26 :			- (a)*(h)*(f));
27 :			}
28 :
29 :	jhr	566	static Diderot_real_t DetM4x4 (Diderot_Mat4x4_t m)
30 :			{
31 :			#define M(i) m[(i)>>2].r[(i)&3]
32 :			return (M( 0) * DetM3x3(M( 5), M( 6), M( 7),
33 :			M( 9), M(10), M(11),
34 :			M(13), M(14), M(15))
35 :			- M( 1) * DetM3x3(M( 4), M( 6), M( 7),
36 :			M( 8), M(10), M(11),
37 :			M(12), M(14), M(15))
38 :			+ M( 2) * DetM3x3(M( 4), M( 5), M( 7),
39 :			M( 8), M( 9), M(11),
40 :			M(12), M(13), M(15))
41 :			- M( 3) * DetM3x3(M( 4), M( 5), M( 6),
42 :			M( 8), M( 9), M(10),
43 :			M(12), M(13), M(14)));
44 :			#undef M
45 :	jhr	439
46 :	jhr	443	}
47 :
48 :			/*! \brief compute the inverse of \arg m, storing the result in \arg i.
49 :			* \param m the matrix to invert
50 :			* \param i the inverted matrix
51 :			*/
52 :	jhr	566	static void InvertM4x4 (Diderot_Mat4x4_t i, Diderot_Mat4x4_t m)
53 :	jhr	439	{
54 :			Diderot_real_t scale = 1.0 / DetM4x4(m);
55 :
56 :	jhr	566	#define M(i) m[(i)>>2].r[(i)&3]
57 :			#define I(j) i[(j)>>2].r[(j)&3]
58 :			I( 0) = DetM3x3(M( 5),M( 6),M( 7),
59 :			M( 9),M(10),M(11),
60 :			M(13),M(14),M(15)) * scale;
61 :	jhr	439
62 :	jhr	566	I( 1) = -DetM3x3(M( 1),M( 2),M( 3),
63 :			M( 9),M(10),M(11),
64 :			M(13),M(14),M(15)) * scale;
65 :	jhr	439
66 :	jhr	566	I( 2) = DetM3x3(M( 1),M( 2),M( 3),
67 :			M( 5),M( 6),M( 7),
68 :			M(13),M(14),M(15)) * scale;
69 :	jhr	439
70 :	jhr	566	I( 3) = -DetM3x3(M( 1),M( 2),M( 3),
71 :			M( 5),M( 6),M( 7),
72 :			M( 9),M(10),M(11)) * scale;
73 :	jhr	439
74 :	jhr	566	I( 4) = -DetM3x3(M( 4),M( 6),M( 7),
75 :			M( 8),M(10),M(11),
76 :			M(12),M(14),M(15)) * scale;
77 :	jhr	439
78 :	jhr	566	I( 5) = DetM3x3(M( 0),M( 2),M( 3),
79 :			M( 8),M(10),M(11),
80 :			M(12),M(14),M(15)) * scale;
81 :	jhr	439
82 :	jhr	566	I( 6) = -DetM3x3(M( 0),M( 2),M( 3),
83 :			M( 4),M( 6),M( 7),
84 :			M(12),M(14),M(15)) * scale;
85 :	jhr	439
86 :	jhr	566	I( 7) = DetM3x3(M( 0),M( 2),M( 3),
87 :			M( 4),M( 6),M( 7),
88 :			M( 8),M(10),M(11)) * scale;
89 :	jhr	439
90 :	jhr	566	I( 8) = DetM3x3(M( 4),M( 5),M( 7),
91 :			M( 8),M( 9),M(11),
92 :			M(12),M(13),M(15)) * scale;
93 :	jhr	439
94 :	jhr	566	I( 9) = -DetM3x3(M( 0),M( 1),M( 3),
95 :			M( 8),M( 9),M(11),
96 :			M(12),M(13),M(15)) * scale;
97 :	jhr	439
98 :	jhr	566	I(10) = DetM3x3(M( 0),M( 1),M( 3),
99 :			M( 4),M( 5),M( 7),
100 :			M(12),M(13),M(15)) * scale;
101 :	jhr	439
102 :	jhr	566	I(11) = -DetM3x3(M( 0),M( 1),M( 3),
103 :			M( 4),M( 5),M( 7),
104 :			M( 8),M( 9),M(11)) * scale;
105 :	jhr	439
106 :	jhr	566	I(12) = -DetM3x3(M( 4),M( 5),M( 6),
107 :			M( 8),M( 9),M(10),
108 :			M(12),M(13),M(14)) * scale;
109 :	jhr	439
110 :	jhr	566	I(13) = DetM3x3(M( 0),M( 1),M( 2),
111 :			M( 8),M( 9),M(10),
112 :			M(12),M(13),M(14)) * scale;
113 :	jhr	439
114 :	jhr	566	I(14) = -DetM3x3(M( 0),M( 1),M( 2),
115 :			M( 4),M( 5),M( 6),
116 :			M(12),M(13),M(14)) * scale;
117 :	jhr	439
118 :	jhr	566	I(15) = DetM3x3(M( 0),M( 1),M( 2),
119 :			M( 4),M( 5),M( 6),
120 :			M( 8),M( 9),M(10)) * scale;
121 :			#undef M
122 :			#undef I
123 :	jhr	439	}
124 :	jhr	443
125 :			/*! \brief compute the inverse transpose of \arg m, storing the result in \arg i.
126 :			* \param m the matrix to invert
127 :			* \param i the inverted matrix
128 :			*/
129 :	jhr	566	static void InvertTransposeM4x4 (Diderot_Mat4x4_t i, Diderot_Mat4x4_t m)
130 :	jhr	443	{
131 :			Diderot_real_t scale = 1.0 / DetM4x4(m);
132 :
133 :	jhr	566	#define M(i) m[(i)>>2].r[(i)&3]
134 :			#define I(j) i[(j)>>2].r[(j)&3]
135 :			I( 0) = DetM3x3(M( 5),M( 6),M( 7),
136 :			M( 9),M(10),M(11),
137 :			M(13),M(14),M(15)) * scale;
138 :	jhr	443
139 :	jhr	566	I( 4) = -DetM3x3(M( 1),M( 2),M( 3),
140 :			M( 9),M(10),M(11),
141 :			M(13),M(14),M(15)) * scale;
142 :	jhr	443
143 :	jhr	566	I( 8) = DetM3x3(M( 1),M( 2),M( 3),
144 :			M( 5),M( 6),M( 7),
145 :			M(13),M(14),M(15)) * scale;
146 :	jhr	443
147 :	jhr	566	I(12) = -DetM3x3(M( 1),M( 2),M( 3),
148 :			M( 5),M( 6),M( 7),
149 :			M( 9),M(10),M(11)) * scale;
150 :	jhr	443
151 :	jhr	566	I( 1) = -DetM3x3(M( 4),M( 6),M( 7),
152 :			M( 8),M(10),M(11),
153 :			M(12),M(14),M(15)) * scale;
154 :	jhr	443
155 :	jhr	566	I( 5) = DetM3x3(M( 0),M( 2),M( 3),
156 :			M( 8),M(10),M(11),
157 :			M(12),M(14),M(15)) * scale;
158 :	jhr	443
159 :	jhr	566	I( 9) = -DetM3x3(M( 0),M( 2),M( 3),
160 :			M( 4),M( 6),M( 7),
161 :			M(12),M(14),M(15)) * scale;
162 :	jhr	443
163 :	jhr	566	I( 7) = DetM3x3(M( 0),M( 2),M( 3),
164 :			M( 4),M( 6),M( 7),
165 :			M( 8),M(10),M(11)) * scale;
166 :	jhr	443
167 :	jhr	566	I( 2) = DetM3x3(M( 4),M( 5),M( 7),
168 :			M( 8),M( 9),M(11),
169 :			M(12),M(13),M(15)) * scale;
170 :	jhr	443
171 :	jhr	566	I( 6) = -DetM3x3(M( 0),M( 1),M( 3),
172 :			M( 8),M( 9),M(11),
173 :			M(12),M(13),M(15)) * scale;
174 :	jhr	443
175 :	jhr	566	I(10) = DetM3x3(M( 0),M( 1),M( 3),
176 :			M( 4),M( 5),M( 7),
177 :			M(12),M(13),M(15)) * scale;
178 :	jhr	443
179 :	jhr	566	I(14) = -DetM3x3(M( 0),M( 1),M( 3),
180 :			M( 4),M( 5),M( 7),
181 :			M( 8),M( 9),M(11)) * scale;
182 :	jhr	443
183 :	jhr	566	I( 3) = -DetM3x3(M( 4),M( 5),M( 6),
184 :			M( 8),M( 9),M(10),
185 :			M(12),M(13),M(14)) * scale;
186 :	jhr	443
187 :	jhr	566	I( 7) = DetM3x3(M( 0),M( 1),M( 2),
188 :			M( 8),M( 9),M(10),
189 :			M(12),M(13),M(14)) * scale;
190 :	jhr	443
191 :	jhr	566	I(11) = -DetM3x3(M( 0),M( 1),M( 2),
192 :			M( 4),M( 5),M( 6),
193 :			M(12),M(13),M(14)) * scale;
194 :	jhr	443
195 :	jhr	566	I(15) = DetM3x3(M( 0),M( 1),M( 2),
196 :			M( 4),M( 5),M( 6),
197 :			M( 8),M( 9),M(10)) * scale;
198 :	jhr	443
199 :	jhr	566	#undef M
200 :			#undef I
201 :	jhr	443	}
202 :	jhr	566
203 :	jhr	568	static void LoadTransformMatrix (Nrrd *nin, Diderot_Mat4x4_t m)
204 :			{
205 :			#define M(i) m[(i)>>2].r[(i)&3]
206 :			int size = nin->spaceDim;
207 :
208 :			// Image axis Scaling and Rotation
209 :			for (int i = 0; i < size; i++) {
210 :			NrrdAxisInfo axisInfo = nin->axis[i];
211 :			for (int j = 0; j < size; j++) {
212 :			int k = (size + 1) * j + i;
213 :			M(k) = axisInfo.spaceDirection[j];
214 :			}
215 :
216 :			//Image Location
217 :			M((i * (size + 1)) + size) = nin->spaceOrigin[i];
218 :
219 :			//Bottom row of the Transform Matrix
220 :			M(((size + 1) * (size)) + i) = 0;
221 :			}
222 :			M(((size + 1) * (size)) + size) = 1;
223 :
224 :			#undef M
225 :			}
226 :
227 :			static Nrrd *LoadNrrdFile (const char *filename)
228 :			{
229 :			/* create a nrrd; at this point this is just an empty container */
230 :			Nrrd *nin = nrrdNew();
231 :
232 :			/* read in the nrrd from file */
233 :			if (nrrdLoad(nin, filename, NULL)) {
234 :			char *err = biffGetDone(NRRD);
235 :			fprintf (stderr, "Diderot: trouble reading \"%s\":\n%s", filename, err);
236 :			free (err);
237 :			return NULL;
238 :			}
239 :
240 :			return nin;
241 :
242 :			}
243 :
244 :	jhr	566	/* load image data from Nrrd files */
245 :			//Status_t Diderot_LoadImage1D (Diderot_string_t name, Diderot_image1D_t *img);
246 :			//Status_t Diderot_LoadImage2D (Diderot_string_t name, Diderot_image2D_t *img);
247 :
248 :			Status_t Diderot_LoadImage3D (Diderot_string_t name, Diderot_image3D_t **imgOut)
249 :			{
250 :	jhr	568	Nrrd *nin = LoadNrrdFile (name);
251 :	jhr	566	Diderot_image3D_t *img = (Diderot_image3D_t *)malloc(sizeof(Diderot_image3D_t));
252 :
253 :			img->dim = 3;
254 :			img->size[0] = nin->axis[0].size;
255 :			img->size[1] = nin->axis[1].size;
256 :			img->size[2] = nin->axis[2].size;
257 :			img->data = nin->data;
258 :
259 :	jhr	568	LoadTransformMatrix (nin, img->m);
260 :	jhr	566	InvertM4x4 (img->m, img->mInv);
261 :			InvertTransposeM4x4 (img->m, img->mInvT);
262 :
263 :			*imgOut = img;
264 :			return DIDEROT_OK;
265 :			}
266 :
267 :			/* functions to get input-parameter values */
268 :			Status_t Diderot_InputString (const char *name, const char **v, bool hasDflt)
269 :			{
270 :	jhr	571	return DIDEROT_FAIL;
271 :	jhr	566	}
272 :
273 :			Status_t Diderot_Inputf (const char *name, float *v, bool hasDflt)
274 :			{
275 :	jhr	571	float f;
276 :
277 :			while (true) {
278 :			if (hasDflt)
279 :			printf("Enter value for %s (default %f): ", name, *v);
280 :			else
281 :			printf("Enter value for %s: ", name);
282 :			fflush (stdout);
283 :
284 :			int n = scanf("%f\n", &f);
285 :			printf("n = %d\n", n);
286 :			if (n == EOF)
287 :			return DIDEROT_FAIL;
288 :			else if (n == 1) {
289 :			*v = f;
290 :			return DIDEROT_OK;;
291 :			}
292 :			else if (hasDflt)
293 :			return DIDEROT_OK;;
294 :			}
295 :
296 :	jhr	566	}
297 :
298 :			Status_t Diderot_InputVec3f (const char *name, vec3f_t *v, bool hasDflt)
299 :			{
300 :	jhr	571	return DIDEROT_FAIL;
301 :	jhr	566	}

---

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