Annotation of /branches/pure-cfg/src/lib/common/image.c

Parent Directory | Revision Log

---

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

1 :	jhr	1093	/*! \file image.c
2 :			*
3 :			* \author John Reppy
4 :	jhr	1213	*
5 :			* \brief code to support the loading of image data from Nrrd files.
6 :	jhr	1093	*/
7 :
8 :			/*
9 :	jhr	1213	* COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
10 :	jhr	1093	* All rights reserved.
11 :			*/
12 :
13 :			#include "Diderot/diderot.h"
14 :			#include <teem/nrrd.h>
15 :
16 :	jhr	1213	/* FIXME: eventally we should change the naming conventions in the header files to be generic */
17 :			#if defined(DIDEROT_SINGLE_PRECISION)
18 :	jhr	1268	#define vec2(a,b) vec2f(a,b)
19 :			#define vec3(a,b,c) vec3f(a,b,c)
20 :			#define mulMat2x2Vec2(m,v) mulMat2x2Vec2f(m,v)
21 :			#define mulMat3x3Vec3(m,v) mulMat3x3Vec3f(m,v)
22 :	jhr	1213	#else
23 :	jhr	1268	#define vec2(a,b) vec2d(a,b)
24 :			#define vec3(a,b,c) vec3d(a,b,c)
25 :			#define mulMat2x2Vec2(m,v) mulMat2x2Vec2d(m,v)
26 :			#define mulMat3x3Vec3(m,v) mulMat3x3Vec3d(m,v)
27 :	jhr	1213	#endif
28 :	jhr	1093
29 :	jhr	1213	// we use double precision representations of the transforms here to
30 :			// improve robustness, even when the results are stored as single-precision
31 :			typedef double Matrix2x2_t[4];
32 :			typedef double Matrix3x3_t[9];
33 :			typedef double Matrix4x4_t[16];
34 :
35 :			typedef double Vec2_t[2];
36 :			typedef double Vec3_t[3];
37 :
38 :			static void PrintMat2x2 (const char *title, Matrix2x2_t m)
39 :	jhr	1093	{
40 :	jhr	1268	#define M(i) m[i]
41 :	jhr	1093	printf ("%s @ %p:\n", title, (void *)m);
42 :			for (int i = 0; i < 2; i++) {
43 :	jhr	1268	int j = 2*i;
44 :			printf (" [ %6lf %6lf ]\n", M(j+0), M(j+1));
45 :	jhr	1093	}
46 :			#undef M
47 :
48 :			}
49 :
50 :	jhr	1213	static void PrintMat3x3 (const char *title, Matrix3x3_t m)
51 :	jhr	1093	{
52 :	jhr	1268	#define M(i) m[i]
53 :	jhr	1093	printf ("%s @ %p:\n", title, (void *)m);
54 :			for (int i = 0; i < 3; i++) {
55 :	jhr	1268	int j = 3*i;
56 :			printf (" [ %6lf %6lf %6lf ]\n", M(j+0), M(j+1), M(j+2));
57 :	jhr	1093	}
58 :			#undef M
59 :
60 :			}
61 :
62 :	jhr	1213	static void PrintMat4x4 (const char *title, Matrix4x4_t m)
63 :	jhr	1093	{
64 :	jhr	1268	#define M(i) m[i]
65 :	jhr	1093	printf ("%s @ %p:\n", title, (void *)m);
66 :			for (int i = 0; i < 4; i++) {
67 :	jhr	1268	int j = 4*i;
68 :			printf (" [ %6lf %6lf %6lf %6lf ]\n", M(j+0), M(j+1), M(j+2), M(j+3));
69 :	jhr	1093	}
70 :			#undef M
71 :
72 :			}
73 :
74 :			/* Loading transfomation information from Nrrd files */
75 :
76 :	jhr	1213	static void LoadTransform1D (Nrrd *nin, double *s, double *t)
77 :	jhr	1093	{
78 :			// compute the offset to the first space axis
79 :			int base = nin->dim - nin->spaceDim;
80 :
81 :			// Image axis Scaling and Rotation
82 :			*s = nin->axis[base+0].spaceDirection[0];
83 :
84 :			// Image location
85 :			*t = nin->spaceOrigin[0];
86 :			}
87 :
88 :	jhr	1213	static void LoadTransform2D (Nrrd *nin, Matrix2x2_t m, Vec2_t t)
89 :	jhr	1093	{
90 :	jhr	1268	#define M(i,j) m[(i)*2+(j)]
91 :	jhr	1093
92 :			// compute the offset to the first space axis
93 :			int base = nin->dim - nin->spaceDim;
94 :
95 :			// Image axis Scaling and Rotation
96 :			M(0,0) = nin->axis[base+0].spaceDirection[0];
97 :			M(1,0) = nin->axis[base+0].spaceDirection[1];
98 :			M(0,1) = nin->axis[base+1].spaceDirection[0];
99 :			M(1,1) = nin->axis[base+1].spaceDirection[1];
100 :
101 :			// Image location
102 :	jhr	1213	t[0] = nin->spaceOrigin[0];
103 :			t[1] = nin->spaceOrigin[1];
104 :	jhr	1093
105 :			#undef M
106 :			}
107 :
108 :	jhr	1213	static void LoadTransform3D (Nrrd *nin, Matrix3x3_t m, Vec3_t t)
109 :	jhr	1093	{
110 :	jhr	1268	#define M(i,j) m[3*(i)+(j)]
111 :	jhr	1093
112 :			// compute the offset to the first space axis
113 :			int base = nin->dim - nin->spaceDim;
114 :
115 :			// Image axis Scaling and Rotation
116 :			M(0,0) = nin->axis[base+0].spaceDirection[0];
117 :			M(1,0) = nin->axis[base+0].spaceDirection[1];
118 :			M(2,0) = nin->axis[base+0].spaceDirection[2];
119 :			M(0,1) = nin->axis[base+1].spaceDirection[0];
120 :			M(1,1) = nin->axis[base+1].spaceDirection[1];
121 :			M(2,1) = nin->axis[base+1].spaceDirection[2];
122 :			M(0,2) = nin->axis[base+2].spaceDirection[0];
123 :			M(1,2) = nin->axis[base+2].spaceDirection[1];
124 :			M(2,2) = nin->axis[base+2].spaceDirection[2];
125 :
126 :			// Image location
127 :	jhr	1213	t[0] = nin->spaceOrigin[0];
128 :			t[1] = nin->spaceOrigin[1];
129 :			t[2] = nin->spaceOrigin[2];
130 :	jhr	1093
131 :			#undef M
132 :			}
133 :
134 :
135 :			/* Transformation matrix operations */
136 :
137 :	jhr	1213	STATIC_INLINE double Det2 (
138 :			double a, double b,
139 :			double c, double d)
140 :	jhr	1093	{
141 :			return (a*d - b*c);
142 :			}
143 :
144 :	jhr	1213	STATIC_INLINE double Det3 (
145 :			double a, double b, double c,
146 :			double d, double e, double f,
147 :			double g, double h, double i)
148 :	jhr	1093	{
149 :			return (a*e*i
150 :	jhr	1268	+ d*h*c
151 :			+ g*b*f
152 :			- g*e*c
153 :			- d*b*i
154 :			- a*h*f);
155 :	jhr	1093	}
156 :
157 :	jhr	1213	static double DetM3x3 (Matrix3x3_t m)
158 :	jhr	1093	{
159 :	jhr	1268	#define M(i) m[i]
160 :	jhr	1093	return Det3(M(0), M(1), M(2),
161 :	jhr	1268	M(3), M(4), M(5),
162 :			M(6), M(7), M(8));
163 :	jhr	1093	#undef M
164 :			}
165 :
166 :	jhr	1213	static double DetM4x4 (Matrix4x4_t m)
167 :	jhr	1093	{
168 :	jhr	1268	#define M(i) m[i]
169 :			return (M( 0) * Det3(M( 5), M( 6), M( 7),
170 :			M( 9), M(10), M(11),
171 :			M(13), M(14), M(15))
172 :			- M( 1) * Det3(M( 4), M( 6), M( 7),
173 :			M( 8), M(10), M(11),
174 :			M(12), M(14), M(15))
175 :			+ M( 2) * Det3(M( 4), M( 5), M( 7),
176 :			M( 8), M( 9), M(11),
177 :			M(12), M(13), M(15))
178 :			- M( 3) * Det3(M( 4), M( 5), M( 6),
179 :			M( 8), M( 9), M(10),
180 :			M(12), M(13), M(14)));
181 :	jhr	1093	#undef M
182 :
183 :			}
184 :
185 :			/*! \brief compute the inverse of \arg m, storing the result in \arg i.
186 :			* \param m the matrix to invert
187 :			* \param i the inverted matrix
188 :			*/
189 :	jhr	1213	static void InvertM2x2 (Matrix2x2_t i, Matrix2x2_t m)
190 :	jhr	1093	{
191 :	jhr	1268	#define M(i) m[i]
192 :			#define I(j) i[j]
193 :			double scale = 1.0 / Det2(M(0), M(1), M(2), M(3));
194 :	jhr	1093
195 :	jhr	1268	I(0) = scale * M(3);
196 :			I(1) = -scale * M(1);
197 :			I(2) = -scale * M(2);
198 :			I(3) = scale * M(0);
199 :	jhr	1093	#undef M
200 :			#undef I
201 :			}
202 :
203 :			/*! \brief compute the inverse of \arg m, storing the result in \arg i.
204 :			* \param m the matrix to invert
205 :			* \param i the inverted matrix
206 :			*/
207 :	jhr	1213	static void InvertM3x3 (Matrix3x3_t i, Matrix3x3_t m)
208 :	jhr	1093	{
209 :	jhr	1268	double scale = 1.0 / DetM3x3(m);
210 :	jhr	1093
211 :	jhr	1268	#define M(i) m[i]
212 :			#define I(j) i[j]
213 :			I(0) = scale * Det2 (M(4), M(5), M(7), M(8));
214 :			I(1) = scale * Det2 (M(2), M(1), M(8), M(7));
215 :			I(2) = scale * Det2 (M(1), M(2), M(4), M(5));
216 :			I(3) = scale * Det2 (M(5), M(3), M(8), M(6));
217 :			I(4) = scale * Det2 (M(0), M(2), M(6), M(8));
218 :			I(5) = scale * Det2 (M(2), M(0), M(5), M(3));
219 :			I(6) = scale * Det2 (M(3), M(4), M(6), M(7));
220 :			I(7) = scale * Det2 (M(1), M(0), M(7), M(6));
221 :			I(8) = scale * Det2 (M(0), M(1), M(3), M(4));
222 :	jhr	1093	#undef M
223 :			#undef I
224 :			}
225 :
226 :			/*! \brief compute the inverse of \arg m, storing the result in \arg i.
227 :			* \param m the matrix to invert
228 :			* \param i the inverted matrix
229 :			*/
230 :	jhr	1213	static void InvertM4x4 (Matrix4x4_t i, Matrix4x4_t m)
231 :	jhr	1093	{
232 :	jhr	1268	double scale = 1.0 / DetM4x4(m);
233 :	jhr	1093
234 :	jhr	1268	#define M(i) m[i]
235 :			#define I(j) i[j]
236 :			I( 0) = Det3(M( 5),M( 6),M( 7),
237 :			M( 9),M(10),M(11),
238 :			M(13),M(14),M(15)) * scale;
239 :
240 :			I( 1) = -Det3(M( 1),M( 2),M( 3),
241 :			M( 9),M(10),M(11),
242 :			M(13),M(14),M(15)) * scale;
243 :
244 :			I( 2) = Det3(M( 1),M( 2),M( 3),
245 :			M( 5),M( 6),M( 7),
246 :			M(13),M(14),M(15)) * scale;
247 :
248 :			I( 3) = -Det3(M( 1),M( 2),M( 3),
249 :			M( 5),M( 6),M( 7),
250 :			M( 9),M(10),M(11)) * scale;
251 :
252 :			I( 4) = -Det3(M( 4),M( 6),M( 7),
253 :			M( 8),M(10),M(11),
254 :			M(12),M(14),M(15)) * scale;
255 :
256 :			I( 5) = Det3(M( 0),M( 2),M( 3),
257 :			M( 8),M(10),M(11),
258 :			M(12),M(14),M(15)) * scale;
259 :
260 :			I( 6) = -Det3(M( 0),M( 2),M( 3),
261 :			M( 4),M( 6),M( 7),
262 :			M(12),M(14),M(15)) * scale;
263 :
264 :			I( 7) = Det3(M( 0),M( 2),M( 3),
265 :			M( 4),M( 6),M( 7),
266 :			M( 8),M(10),M(11)) * scale;
267 :
268 :			I( 8) = Det3(M( 4),M( 5),M( 7),
269 :			M( 8),M( 9),M(11),
270 :			M(12),M(13),M(15)) * scale;
271 :
272 :			I( 9) = -Det3(M( 0),M( 1),M( 3),
273 :			M( 8),M( 9),M(11),
274 :			M(12),M(13),M(15)) * scale;
275 :
276 :			I(10) = Det3(M( 0),M( 1),M( 3),
277 :			M( 4),M( 5),M( 7),
278 :			M(12),M(13),M(15)) * scale;
279 :
280 :			I(11) = -Det3(M( 0),M( 1),M( 3),
281 :			M( 4),M( 5),M( 7),
282 :			M( 8),M( 9),M(11)) * scale;
283 :
284 :			I(12) = -Det3(M( 4),M( 5),M( 6),
285 :			M( 8),M( 9),M(10),
286 :			M(12),M(13),M(14)) * scale;
287 :
288 :			I(13) = Det3(M( 0),M( 1),M( 2),
289 :			M( 8),M( 9),M(10),
290 :			M(12),M(13),M(14)) * scale;
291 :
292 :			I(14) = -Det3(M( 0),M( 1),M( 2),
293 :			M( 4),M( 5),M( 6),
294 :			M(12),M(13),M(14)) * scale;
295 :
296 :			I(15) = Det3(M( 0),M( 1),M( 2),
297 :			M( 4),M( 5),M( 6),
298 :			M( 8),M( 9),M(10)) * scale;
299 :	jhr	1093	#undef M
300 :			#undef I
301 :			}
302 :
303 :			static Nrrd *LoadNrrdFile (const char *filename)
304 :			{
305 :			/* create a nrrd; at this point this is just an empty container */
306 :			Nrrd *nin = nrrdNew();
307 :
308 :			/* read in the nrrd from file */
309 :			if (nrrdLoad(nin, filename, NULL)) {
310 :	jhr	1268	char *err = biffGetDone(NRRD);
311 :			fprintf (stderr, "Diderot: trouble reading \"%s\":\n%s", filename, err);
312 :			free (err);
313 :			return NULL;
314 :	jhr	1093	}
315 :
316 :			return nin;
317 :
318 :			}
319 :
320 :			STATIC_INLINE bool CheckAxis (Nrrd *nin, int i)
321 :			{
322 :			return (nin->axis[i].kind != nrrdKindSpace);
323 :			}
324 :
325 :			/* load image data from Nrrd files */
326 :			Status_t Diderot_LoadImage1D (Diderot_string_t name, Diderot_image1D_t **imgOut)
327 :			{
328 :			Nrrd *nin = LoadNrrdFile (name);
329 :			Diderot_image1D_t *img = (Diderot_image1D_t *)malloc(sizeof(Diderot_image1D_t));
330 :
331 :			// compute the offset to the first space axis
332 :			int base = nin->dim - nin->spaceDim;
333 :
334 :	jhr	1357	// FIXME: we should also be checking the shape of the voxels and the type of the samples.
335 :	jhr	1093	// check that the image format is what we expect
336 :			if (nin->spaceDim != 1) {
337 :	jhr	1268	fprintf(stderr, "nrrd file \"%s\" has unexpected dimension %d, expected 1\n",
338 :			name, nin->spaceDim);
339 :			return DIDEROT_FAIL;
340 :	jhr	1093	}
341 :			else if (CheckAxis(nin, base)) {
342 :	jhr	1268	fprintf(stderr,
343 :			"nrrd file \"%s\" has unexpected axis structure: %s\n",
344 :			name,
345 :			airEnumStr(nrrdKind, nin->axis[base].kind));
346 :			return DIDEROT_FAIL;
347 :	jhr	1093	}
348 :
349 :			img->dim = 1;
350 :			img->size[0] = nin->axis[base+0].size;
351 :			img->data = nin->data;
352 :	jhr	1357	// FIXME: this code probably only works for scalar images
353 :			img->dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin) * img->size[0];
354 :	jhr	1093
355 :	jhr	1273	//printf("LoadImage \"%s\": space dim = 1, axes = <%d>\n",
356 :			//name, img->size[base+0]);
357 :	jhr	1093
358 :			// from the Nrrd file, we load the scaling and translation
359 :	jhr	1213	double s, t;
360 :	jhr	1093	LoadTransform1D (nin, &s, &t);
361 :
362 :	jhr	1213	img->s = (Diderot_real_t)(1.0 / s);
363 :			img->t = (Diderot_real_t)(-t / s);
364 :	jhr	1093
365 :			*imgOut = img;
366 :			return DIDEROT_OK;
367 :			}
368 :
369 :			Status_t Diderot_LoadImage2D (Diderot_string_t name, Diderot_image2D_t **imgOut)
370 :			{
371 :			Nrrd *nin = LoadNrrdFile (name);
372 :			Diderot_image2D_t *img = (Diderot_image2D_t *)malloc(sizeof(Diderot_image2D_t));
373 :
374 :			// compute the offset to the first space axis
375 :			int base = nin->dim - nin->spaceDim;
376 :
377 :			// check that the image format is what we expect
378 :			if (nin->spaceDim != 2) {
379 :	jhr	1268	fprintf(stderr, "nrrd file \"%s\" has unexpected dimension %d, expected 2\n",
380 :			name, nin->spaceDim);
381 :			return DIDEROT_FAIL;
382 :	jhr	1093	}
383 :			else if (CheckAxis(nin, base) || CheckAxis(nin, base+1)) {
384 :	jhr	1268	fprintf(stderr,
385 :			"nrrd file \"%s\" has unexpected axis structure: %s %s\n",
386 :			name,
387 :			airEnumStr(nrrdKind, nin->axis[base].kind),
388 :			airEnumStr(nrrdKind, nin->axis[base+1].kind));
389 :			return DIDEROT_FAIL;
390 :	jhr	1093	}
391 :
392 :			img->dim = 2;
393 :			img->size[0] = nin->axis[base+0].size;
394 :			img->size[1] = nin->axis[base+1].size;
395 :			img->data = nin->data;
396 :	jhr	1357	// FIXME: this code probably only works for scalar images
397 :			img->dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin) * img->size[0] * img->size[1];
398 :	jhr	1093
399 :	jhr	1273	//printf("LoadImage \"%s\": space dim = 2, axes = <%d, %d>\n",
400 :			//name, img->size[0], img->size[1]);
401 :	jhr	1093
402 :			// from the Nrrd file, we load the affine image-to-world transform matrix
403 :	jhr	1268	Matrix2x2_t m; // rotation and scaling
404 :			Vec2_t t; // translation part of the transform
405 :	jhr	1213	LoadTransform2D (nin, m, t);
406 :	jhr	1093
407 :	jhr	1273	//PrintMat2x2("m", m);
408 :			//printf("t = <%f, %f>\n", t[0], t[1]);
409 :	jhr	1093
410 :			// compute inverse of m, which is the transform from the world basis to the image basis
411 :	jhr	1213	Matrix2x2_t mInv;
412 :			InvertM2x2 (mInv, m);
413 :	jhr	1273	//PrintMat2x2("w2i", mInv);
414 :	jhr	1093
415 :	jhr	1213	// copy results into the image data structure
416 :			for (int i = 0; i < 2; i++) {
417 :	jhr	1268	for (int j = 0; j < 2; j++) {
418 :			img->w2i[i].r[j] = mInv[2*i+j];
419 :			}
420 :	jhr	1213	}
421 :	jhr	1093	transpose2x2f (img->w2iT, img->w2i);
422 :	jhr	1213	Diderot_vec2_t tVec = vec2(t[0], t[1]);
423 :			img->tVec = -mulMat2x2Vec2(img->w2i, tVec); // inv([M t]) = [inv(M) -inv(M)t]
424 :	jhr	1273	//printf("tVec = <%f, %f>\n", ((union2f_t)img->tVec).r[0], ((union2f_t)img->tVec).r[1]);
425 :	jhr	1093
426 :			*imgOut = img;
427 :			return DIDEROT_OK;
428 :			}
429 :
430 :			Status_t Diderot_LoadImage3D (Diderot_string_t name, Diderot_image3D_t **imgOut)
431 :			{
432 :			Nrrd *nin = LoadNrrdFile (name);
433 :
434 :			// compute the offset to the first space axis
435 :			int base = nin->dim - nin->spaceDim;
436 :
437 :			// check that the image format is what we expect
438 :			if (nin->spaceDim != 3) {
439 :	jhr	1268	fprintf(stderr, "nrrd file \"%s\" has unexpected dimension %d, expected 3\n",
440 :			name, nin->spaceDim);
441 :			return DIDEROT_FAIL;
442 :	jhr	1093	}
443 :			else if (CheckAxis(nin, base) || CheckAxis(nin, base+1) || CheckAxis(nin, base+2)) {
444 :	jhr	1268	fprintf(stderr,
445 :			"nrrd file \"%s\" has unexpected axis structure: %s %s %s\n",
446 :			name,
447 :			airEnumStr(nrrdKind, nin->axis[base].kind),
448 :			airEnumStr(nrrdKind, nin->axis[base+1].kind),
449 :			airEnumStr(nrrdKind, nin->axis[base+2].kind));
450 :			return DIDEROT_FAIL;
451 :	jhr	1093	}
452 :
453 :			Diderot_image3D_t *img = (Diderot_image3D_t *)malloc(sizeof(Diderot_image3D_t));
454 :
455 :			img->dim = 3;
456 :			img->size[0] = nin->axis[base+0].size;
457 :			img->size[1] = nin->axis[base+1].size;
458 :			img->size[2] = nin->axis[base+2].size;
459 :			img->data = nin->data;
460 :	jhr	1357	// FIXME: this code probably only works for scalar images
461 :			img->dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin) * img->size[0] * img->size[1] * img->size[2];
462 :	jhr	1093
463 :	jhr	1273	//printf("LoadImage \"%s\": space dim = 3, axes = <%d, %d, %d>\n",
464 :			//name, img->size[0], img->size[1], img->size[2]);
465 :	jhr	1093
466 :			// from the Nrrd file, we load the affine image-to-world transform matrix
467 :	jhr	1268	Matrix3x3_t m; // rotation and scaling
468 :			Vec3_t t; // translation part of the transform
469 :	jhr	1213	LoadTransform3D (nin, m, t);
470 :	jhr	1093
471 :	jhr	1273	//PrintMat3x3("m", m);
472 :			//printf("t = <%f, %f, %f>\n", t[0], t[1], t[2]);
473 :	jhr	1093
474 :			// compute inverse of m, which is the transform from the world basis to the image basis
475 :	jhr	1213	Matrix3x3_t mInv;
476 :			InvertM3x3 (mInv, m);
477 :	jhr	1273	//PrintMat3x3("w2i", mInv);
478 :	jhr	1093
479 :	jhr	1213	// copy results into the image data structure
480 :			for (int i = 0; i < 3; i++) {
481 :	jhr	1268	for (int j = 0; j < 3; j++) {
482 :			img->w2i[i].r[j] = mInv[3*i+j];
483 :			}
484 :	jhr	1213	}
485 :	jhr	1093	transpose3x3f (img->w2iT, img->w2i);
486 :	jhr	1213	Diderot_vec3_t tVec = vec3(t[0], t[1], t[2]);
487 :			img->tVec = -mulMat3x3Vec3(img->w2i, tVec); // inv([M t]) = [inv(M) -inv(M)t]
488 :	jhr	1273	//printf("tVec = <%f, %f, %f>\n", ((Diderot_union3_t)img->tVec).r[0], ((Diderot_union3_t)img->tVec).r[1], ((Diderot_union3_t)img->tVec).r[2]);
489 :	jhr	1093
490 :			*imgOut = img;
491 :			return DIDEROT_OK;
492 :			}

---

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