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

Parent Directory | Revision Log

---

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

---

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