Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /branches/pure-cfg/src/lib/common/image.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1364 - (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 1364 img->dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin);
353 : jhr 1093
354 : jhr 1273 //printf("LoadImage \"%s\": space dim = 1, axes = <%d>\n",
355 :     //name, img->size[base+0]);
356 : jhr 1093
357 :     // from the Nrrd file, we load the scaling and translation
358 : jhr 1213 double s, t;
359 : jhr 1093 LoadTransform1D (nin, &s, &t);
360 :    
361 : jhr 1213 img->s = (Diderot_real_t)(1.0 / s);
362 :     img->t = (Diderot_real_t)(-t / s);
363 : jhr 1093
364 :     *imgOut = img;
365 :     return DIDEROT_OK;
366 :     }
367 :    
368 :     Status_t Diderot_LoadImage2D (Diderot_string_t name, Diderot_image2D_t **imgOut)
369 :     {
370 :     Nrrd *nin = LoadNrrdFile (name);
371 :     Diderot_image2D_t *img = (Diderot_image2D_t *)malloc(sizeof(Diderot_image2D_t));
372 :    
373 :     // compute the offset to the first space axis
374 :     int base = nin->dim - nin->spaceDim;
375 :    
376 :     // check that the image format is what we expect
377 :     if (nin->spaceDim != 2) {
378 : jhr 1268 fprintf(stderr, "nrrd file \"%s\" has unexpected dimension %d, expected 2\n",
379 :     name, nin->spaceDim);
380 :     return DIDEROT_FAIL;
381 : jhr 1093 }
382 :     else if (CheckAxis(nin, base) || CheckAxis(nin, base+1)) {
383 : jhr 1268 fprintf(stderr,
384 :     "nrrd file \"%s\" has unexpected axis structure: %s %s\n",
385 :     name,
386 :     airEnumStr(nrrdKind, nin->axis[base].kind),
387 :     airEnumStr(nrrdKind, nin->axis[base+1].kind));
388 :     return DIDEROT_FAIL;
389 : jhr 1093 }
390 :    
391 :     img->dim = 2;
392 :     img->size[0] = nin->axis[base+0].size;
393 :     img->size[1] = nin->axis[base+1].size;
394 :     img->data = nin->data;
395 : jhr 1364 img->dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin);
396 : jhr 1093
397 : jhr 1273 //printf("LoadImage \"%s\": space dim = 2, axes = <%d, %d>\n",
398 :     //name, img->size[0], img->size[1]);
399 : jhr 1093
400 :     // from the Nrrd file, we load the affine image-to-world transform matrix
401 : jhr 1268 Matrix2x2_t m; // rotation and scaling
402 :     Vec2_t t; // translation part of the transform
403 : jhr 1213 LoadTransform2D (nin, m, t);
404 : jhr 1093
405 : jhr 1273 //PrintMat2x2("m", m);
406 :     //printf("t = <%f, %f>\n", t[0], t[1]);
407 : jhr 1093
408 :     // compute inverse of m, which is the transform from the world basis to the image basis
409 : jhr 1213 Matrix2x2_t mInv;
410 :     InvertM2x2 (mInv, m);
411 : jhr 1273 //PrintMat2x2("w2i", mInv);
412 : jhr 1093
413 : jhr 1213 // copy results into the image data structure
414 :     for (int i = 0; i < 2; i++) {
415 : jhr 1268 for (int j = 0; j < 2; j++) {
416 :     img->w2i[i].r[j] = mInv[2*i+j];
417 :     }
418 : jhr 1213 }
419 : jhr 1093 transpose2x2f (img->w2iT, img->w2i);
420 : jhr 1213 Diderot_vec2_t tVec = vec2(t[0], t[1]);
421 :     img->tVec = -mulMat2x2Vec2(img->w2i, tVec); // inv([M t]) = [inv(M) -inv(M)t]
422 : jhr 1273 //printf("tVec = <%f, %f>\n", ((union2f_t)img->tVec).r[0], ((union2f_t)img->tVec).r[1]);
423 : jhr 1093
424 :     *imgOut = img;
425 :     return DIDEROT_OK;
426 :     }
427 :    
428 :     Status_t Diderot_LoadImage3D (Diderot_string_t name, Diderot_image3D_t **imgOut)
429 :     {
430 :     Nrrd *nin = LoadNrrdFile (name);
431 :    
432 :     // compute the offset to the first space axis
433 :     int base = nin->dim - nin->spaceDim;
434 :    
435 :     // check that the image format is what we expect
436 :     if (nin->spaceDim != 3) {
437 : jhr 1268 fprintf(stderr, "nrrd file \"%s\" has unexpected dimension %d, expected 3\n",
438 :     name, nin->spaceDim);
439 :     return DIDEROT_FAIL;
440 : jhr 1093 }
441 :     else if (CheckAxis(nin, base) || CheckAxis(nin, base+1) || CheckAxis(nin, base+2)) {
442 : jhr 1268 fprintf(stderr,
443 :     "nrrd file \"%s\" has unexpected axis structure: %s %s %s\n",
444 :     name,
445 :     airEnumStr(nrrdKind, nin->axis[base].kind),
446 :     airEnumStr(nrrdKind, nin->axis[base+1].kind),
447 :     airEnumStr(nrrdKind, nin->axis[base+2].kind));
448 :     return DIDEROT_FAIL;
449 : jhr 1093 }
450 :    
451 :     Diderot_image3D_t *img = (Diderot_image3D_t *)malloc(sizeof(Diderot_image3D_t));
452 :    
453 :     img->dim = 3;
454 :     img->size[0] = nin->axis[base+0].size;
455 :     img->size[1] = nin->axis[base+1].size;
456 :     img->size[2] = nin->axis[base+2].size;
457 :     img->data = nin->data;
458 : jhr 1364 img->dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin);
459 : jhr 1093
460 : jhr 1273 //printf("LoadImage \"%s\": space dim = 3, axes = <%d, %d, %d>\n",
461 :     //name, img->size[0], img->size[1], img->size[2]);
462 : jhr 1093
463 :     // from the Nrrd file, we load the affine image-to-world transform matrix
464 : jhr 1268 Matrix3x3_t m; // rotation and scaling
465 :     Vec3_t t; // translation part of the transform
466 : jhr 1213 LoadTransform3D (nin, m, t);
467 : jhr 1093
468 : jhr 1273 //PrintMat3x3("m", m);
469 :     //printf("t = <%f, %f, %f>\n", t[0], t[1], t[2]);
470 : jhr 1093
471 :     // compute inverse of m, which is the transform from the world basis to the image basis
472 : jhr 1213 Matrix3x3_t mInv;
473 :     InvertM3x3 (mInv, m);
474 : jhr 1273 //PrintMat3x3("w2i", mInv);
475 : jhr 1093
476 : jhr 1213 // copy results into the image data structure
477 :     for (int i = 0; i < 3; i++) {
478 : jhr 1268 for (int j = 0; j < 3; j++) {
479 :     img->w2i[i].r[j] = mInv[3*i+j];
480 :     }
481 : jhr 1213 }
482 : jhr 1093 transpose3x3f (img->w2iT, img->w2i);
483 : jhr 1213 Diderot_vec3_t tVec = vec3(t[0], t[1], t[2]);
484 :     img->tVec = -mulMat3x3Vec3(img->w2i, tVec); // inv([M t]) = [inv(M) -inv(M)t]
485 : 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]);
486 : jhr 1093
487 :     *imgOut = img;
488 :     return DIDEROT_OK;
489 :     }

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