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

SCM Repository

[diderot] Annotation of /trunk/src/lib/common/image.c
ViewVC logotype

Annotation of /trunk/src/lib/common/image.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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

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