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

SCM Repository

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

Annotation of /branches/pure-cfg/src/lib/diderot.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : jhr 439 /*! \file diderot.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 : jhr 566 #include "Diderot/diderot.h"
12 : jhr 439 #include <teem/nrrd.h>
13 :    
14 :     /* Transformation matrix operations */
15 :    
16 :     STATIC_INLINE Diderot_real_t DetM3x3 (
17 :     Diderot_real_t a, Diderot_real_t b, Diderot_real_t c,
18 :     Diderot_real_t d, Diderot_real_t e, Diderot_real_t f,
19 :     Diderot_real_t g, Diderot_real_t h, Diderot_real_t i)
20 :     {
21 :     return ( (a)*(e)*(i)
22 :     + (d)*(h)*(c)
23 :     + (g)*(b)*(f)
24 :     - (g)*(e)*(c)
25 :     - (d)*(b)*(i)
26 :     - (a)*(h)*(f));
27 :     }
28 :    
29 : jhr 566 static Diderot_real_t DetM4x4 (Diderot_Mat4x4_t m)
30 :     {
31 :     #define M(i) m[(i)>>2].r[(i)&3]
32 :     return (M( 0) * DetM3x3(M( 5), M( 6), M( 7),
33 :     M( 9), M(10), M(11),
34 :     M(13), M(14), M(15))
35 :     - M( 1) * DetM3x3(M( 4), M( 6), M( 7),
36 :     M( 8), M(10), M(11),
37 :     M(12), M(14), M(15))
38 :     + M( 2) * DetM3x3(M( 4), M( 5), M( 7),
39 :     M( 8), M( 9), M(11),
40 :     M(12), M(13), M(15))
41 :     - M( 3) * DetM3x3(M( 4), M( 5), M( 6),
42 :     M( 8), M( 9), M(10),
43 :     M(12), M(13), M(14)));
44 :     #undef M
45 : jhr 439
46 : jhr 443 }
47 :    
48 :     /*! \brief compute the inverse of \arg m, storing the result in \arg i.
49 :     * \param m the matrix to invert
50 :     * \param i the inverted matrix
51 :     */
52 : jhr 566 static void InvertM4x4 (Diderot_Mat4x4_t i, Diderot_Mat4x4_t m)
53 : jhr 439 {
54 :     Diderot_real_t scale = 1.0 / DetM4x4(m);
55 :    
56 : jhr 566 #define M(i) m[(i)>>2].r[(i)&3]
57 :     #define I(j) i[(j)>>2].r[(j)&3]
58 :     I( 0) = DetM3x3(M( 5),M( 6),M( 7),
59 :     M( 9),M(10),M(11),
60 :     M(13),M(14),M(15)) * scale;
61 : jhr 439
62 : jhr 566 I( 1) = -DetM3x3(M( 1),M( 2),M( 3),
63 :     M( 9),M(10),M(11),
64 :     M(13),M(14),M(15)) * scale;
65 : jhr 439
66 : jhr 566 I( 2) = DetM3x3(M( 1),M( 2),M( 3),
67 :     M( 5),M( 6),M( 7),
68 :     M(13),M(14),M(15)) * scale;
69 : jhr 439
70 : jhr 566 I( 3) = -DetM3x3(M( 1),M( 2),M( 3),
71 :     M( 5),M( 6),M( 7),
72 :     M( 9),M(10),M(11)) * scale;
73 : jhr 439
74 : jhr 566 I( 4) = -DetM3x3(M( 4),M( 6),M( 7),
75 :     M( 8),M(10),M(11),
76 :     M(12),M(14),M(15)) * scale;
77 : jhr 439
78 : jhr 566 I( 5) = DetM3x3(M( 0),M( 2),M( 3),
79 :     M( 8),M(10),M(11),
80 :     M(12),M(14),M(15)) * scale;
81 : jhr 439
82 : jhr 566 I( 6) = -DetM3x3(M( 0),M( 2),M( 3),
83 :     M( 4),M( 6),M( 7),
84 :     M(12),M(14),M(15)) * scale;
85 : jhr 439
86 : jhr 566 I( 7) = DetM3x3(M( 0),M( 2),M( 3),
87 :     M( 4),M( 6),M( 7),
88 :     M( 8),M(10),M(11)) * scale;
89 : jhr 439
90 : jhr 566 I( 8) = DetM3x3(M( 4),M( 5),M( 7),
91 :     M( 8),M( 9),M(11),
92 :     M(12),M(13),M(15)) * scale;
93 : jhr 439
94 : jhr 566 I( 9) = -DetM3x3(M( 0),M( 1),M( 3),
95 :     M( 8),M( 9),M(11),
96 :     M(12),M(13),M(15)) * scale;
97 : jhr 439
98 : jhr 566 I(10) = DetM3x3(M( 0),M( 1),M( 3),
99 :     M( 4),M( 5),M( 7),
100 :     M(12),M(13),M(15)) * scale;
101 : jhr 439
102 : jhr 566 I(11) = -DetM3x3(M( 0),M( 1),M( 3),
103 :     M( 4),M( 5),M( 7),
104 :     M( 8),M( 9),M(11)) * scale;
105 : jhr 439
106 : jhr 566 I(12) = -DetM3x3(M( 4),M( 5),M( 6),
107 :     M( 8),M( 9),M(10),
108 :     M(12),M(13),M(14)) * scale;
109 : jhr 439
110 : jhr 566 I(13) = DetM3x3(M( 0),M( 1),M( 2),
111 :     M( 8),M( 9),M(10),
112 :     M(12),M(13),M(14)) * scale;
113 : jhr 439
114 : jhr 566 I(14) = -DetM3x3(M( 0),M( 1),M( 2),
115 :     M( 4),M( 5),M( 6),
116 :     M(12),M(13),M(14)) * scale;
117 : jhr 439
118 : jhr 566 I(15) = DetM3x3(M( 0),M( 1),M( 2),
119 :     M( 4),M( 5),M( 6),
120 :     M( 8),M( 9),M(10)) * scale;
121 :     #undef M
122 :     #undef I
123 : jhr 439 }
124 : jhr 443
125 :     /*! \brief compute the inverse transpose of \arg m, storing the result in \arg i.
126 :     * \param m the matrix to invert
127 :     * \param i the inverted matrix
128 :     */
129 : jhr 566 static void InvertTransposeM4x4 (Diderot_Mat4x4_t i, Diderot_Mat4x4_t m)
130 : jhr 443 {
131 :     Diderot_real_t scale = 1.0 / DetM4x4(m);
132 :    
133 : jhr 566 #define M(i) m[(i)>>2].r[(i)&3]
134 :     #define I(j) i[(j)>>2].r[(j)&3]
135 :     I( 0) = DetM3x3(M( 5),M( 6),M( 7),
136 :     M( 9),M(10),M(11),
137 :     M(13),M(14),M(15)) * scale;
138 : jhr 443
139 : jhr 566 I( 4) = -DetM3x3(M( 1),M( 2),M( 3),
140 :     M( 9),M(10),M(11),
141 :     M(13),M(14),M(15)) * scale;
142 : jhr 443
143 : jhr 566 I( 8) = DetM3x3(M( 1),M( 2),M( 3),
144 :     M( 5),M( 6),M( 7),
145 :     M(13),M(14),M(15)) * scale;
146 : jhr 443
147 : jhr 566 I(12) = -DetM3x3(M( 1),M( 2),M( 3),
148 :     M( 5),M( 6),M( 7),
149 :     M( 9),M(10),M(11)) * scale;
150 : jhr 443
151 : jhr 566 I( 1) = -DetM3x3(M( 4),M( 6),M( 7),
152 :     M( 8),M(10),M(11),
153 :     M(12),M(14),M(15)) * scale;
154 : jhr 443
155 : jhr 566 I( 5) = DetM3x3(M( 0),M( 2),M( 3),
156 :     M( 8),M(10),M(11),
157 :     M(12),M(14),M(15)) * scale;
158 : jhr 443
159 : jhr 566 I( 9) = -DetM3x3(M( 0),M( 2),M( 3),
160 :     M( 4),M( 6),M( 7),
161 :     M(12),M(14),M(15)) * scale;
162 : jhr 443
163 : jhr 566 I( 7) = DetM3x3(M( 0),M( 2),M( 3),
164 :     M( 4),M( 6),M( 7),
165 :     M( 8),M(10),M(11)) * scale;
166 : jhr 443
167 : jhr 566 I( 2) = DetM3x3(M( 4),M( 5),M( 7),
168 :     M( 8),M( 9),M(11),
169 :     M(12),M(13),M(15)) * scale;
170 : jhr 443
171 : jhr 566 I( 6) = -DetM3x3(M( 0),M( 1),M( 3),
172 :     M( 8),M( 9),M(11),
173 :     M(12),M(13),M(15)) * scale;
174 : jhr 443
175 : jhr 566 I(10) = DetM3x3(M( 0),M( 1),M( 3),
176 :     M( 4),M( 5),M( 7),
177 :     M(12),M(13),M(15)) * scale;
178 : jhr 443
179 : jhr 566 I(14) = -DetM3x3(M( 0),M( 1),M( 3),
180 :     M( 4),M( 5),M( 7),
181 :     M( 8),M( 9),M(11)) * scale;
182 : jhr 443
183 : jhr 566 I( 3) = -DetM3x3(M( 4),M( 5),M( 6),
184 :     M( 8),M( 9),M(10),
185 :     M(12),M(13),M(14)) * scale;
186 : jhr 443
187 : jhr 566 I( 7) = DetM3x3(M( 0),M( 1),M( 2),
188 :     M( 8),M( 9),M(10),
189 :     M(12),M(13),M(14)) * scale;
190 : jhr 443
191 : jhr 566 I(11) = -DetM3x3(M( 0),M( 1),M( 2),
192 :     M( 4),M( 5),M( 6),
193 :     M(12),M(13),M(14)) * scale;
194 : jhr 443
195 : jhr 566 I(15) = DetM3x3(M( 0),M( 1),M( 2),
196 :     M( 4),M( 5),M( 6),
197 :     M( 8),M( 9),M(10)) * scale;
198 : jhr 443
199 : jhr 566 #undef M
200 :     #undef I
201 : jhr 443 }
202 : jhr 566
203 : jhr 568 static void LoadTransformMatrix (Nrrd *nin, Diderot_Mat4x4_t m)
204 :     {
205 :     #define M(i) m[(i)>>2].r[(i)&3]
206 :     int size = nin->spaceDim;
207 :    
208 :     // Image axis Scaling and Rotation
209 :     for (int i = 0; i < size; i++) {
210 :     NrrdAxisInfo axisInfo = nin->axis[i];
211 :     for (int j = 0; j < size; j++) {
212 :     int k = (size + 1) * j + i;
213 :     M(k) = axisInfo.spaceDirection[j];
214 :     }
215 :    
216 :     //Image Location
217 :     M((i * (size + 1)) + size) = nin->spaceOrigin[i];
218 :    
219 :     //Bottom row of the Transform Matrix
220 :     M(((size + 1) * (size)) + i) = 0;
221 :     }
222 :     M(((size + 1) * (size)) + size) = 1;
223 :    
224 :     #undef M
225 :     }
226 :    
227 :     static Nrrd *LoadNrrdFile (const char *filename)
228 :     {
229 :     /* create a nrrd; at this point this is just an empty container */
230 :     Nrrd *nin = nrrdNew();
231 :    
232 :     /* read in the nrrd from file */
233 :     if (nrrdLoad(nin, filename, NULL)) {
234 :     char *err = biffGetDone(NRRD);
235 :     fprintf (stderr, "Diderot: trouble reading \"%s\":\n%s", filename, err);
236 :     free (err);
237 :     return NULL;
238 :     }
239 :    
240 :     return nin;
241 :    
242 :     }
243 :    
244 : jhr 566 /* load image data from Nrrd files */
245 :     //Status_t Diderot_LoadImage1D (Diderot_string_t name, Diderot_image1D_t *img);
246 :     //Status_t Diderot_LoadImage2D (Diderot_string_t name, Diderot_image2D_t *img);
247 :    
248 :     Status_t Diderot_LoadImage3D (Diderot_string_t name, Diderot_image3D_t **imgOut)
249 :     {
250 : jhr 568 Nrrd *nin = LoadNrrdFile (name);
251 : jhr 566 Diderot_image3D_t *img = (Diderot_image3D_t *)malloc(sizeof(Diderot_image3D_t));
252 :    
253 :     img->dim = 3;
254 :     img->size[0] = nin->axis[0].size;
255 :     img->size[1] = nin->axis[1].size;
256 :     img->size[2] = nin->axis[2].size;
257 :     img->data = nin->data;
258 :    
259 : jhr 568 LoadTransformMatrix (nin, img->m);
260 : jhr 566 InvertM4x4 (img->m, img->mInv);
261 :     InvertTransposeM4x4 (img->m, img->mInvT);
262 :    
263 :     *imgOut = img;
264 :     return DIDEROT_OK;
265 :     }
266 :    
267 :     /* functions to get input-parameter values */
268 :     Status_t Diderot_InputString (const char *name, const char **v, bool hasDflt)
269 :     {
270 : jhr 571 return DIDEROT_FAIL;
271 : jhr 566 }
272 :    
273 :     Status_t Diderot_Inputf (const char *name, float *v, bool hasDflt)
274 :     {
275 : jhr 572 char buf[256];
276 : jhr 571 float f;
277 :    
278 :     while (true) {
279 :     if (hasDflt)
280 :     printf("Enter value for %s (default %f): ", name, *v);
281 :     else
282 :     printf("Enter value for %s: ", name);
283 :     fflush (stdout);
284 :    
285 : jhr 572 char *p = fgets(buf, sizeof(buf), stdin);
286 :     if (p == NULL)
287 :     return DIDEROT_FAIL; // EOF
288 :     int n = sscanf(buf, "%f\n", &f);
289 :     if (n == 1) {
290 : jhr 571 *v = f;
291 :     return DIDEROT_OK;;
292 :     }
293 :     else if (hasDflt)
294 :     return DIDEROT_OK;;
295 :     }
296 :    
297 : jhr 566 }
298 :    
299 :     Status_t Diderot_InputVec3f (const char *name, vec3f_t *v, bool hasDflt)
300 :     {
301 : jhr 572 char buf[256];
302 :     float f1, f2, f3;
303 :    
304 :     while (true) {
305 :     if (hasDflt)
306 :     printf("Enter value for %s (default %f %f %f): ",
307 :     name, ((union3f_t)*v).r[0], ((union3f_t)*v).r[1], ((union3f_t)*v).r[2]);
308 :     else
309 :     printf("Enter value for %s: ", name);
310 :     fflush (stdout);
311 :    
312 :     char *p = fgets(buf, sizeof(buf), stdin);
313 :     if (p == NULL)
314 :     return DIDEROT_FAIL; // EOF
315 :     int n = sscanf(buf, "%f %f %f\n", &f1, &f2, &f3);
316 :     if (n == 3) {
317 :     *v = vec3f(f1, f2, f3);
318 :     return DIDEROT_OK;;
319 :     }
320 :     else if (hasDflt)
321 :     return DIDEROT_OK;;
322 :     }
323 :    
324 : jhr 566 }

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