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 1093 - (view) (download) (as text)

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

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