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

SCM Repository

[diderot] Annotation of /branches/vis15/src/lib/common/image.cxx
ViewVC logotype

Annotation of /branches/vis15/src/lib/common/image.cxx

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : jhr 3823 /*! \file image.cxx
2 :     *
3 :     * \author John Reppy
4 :     */
5 :    
6 :     /*
7 :     * This code is part of the Diderot Project (http://diderot-language.cs.uchicago.edu)
8 :     *
9 :     * COPYRIGHT (c) 2016 The University of Chicago
10 :     * All rights reserved.
11 :     */
12 :    
13 :     #include <cstdlib>
14 :     #include "image.hxx"
15 :    
16 :     namespace Diderot {
17 :    
18 :     /******************** 1D images ********************/
19 :    
20 :     static void load_transform_1d (const Nrrd *nin, double &s, double &t)
21 :     {
22 :     // compute the offset to the first space axis
23 :     int base = nin->dim - nin->spaceDim;
24 :    
25 :     // Image axis Scaling and Rotation
26 :     s = nin->axis[base+0].spaceDirection[0];
27 :    
28 :     // Image location
29 :     t = nin->spaceOrigin[0];
30 :     }
31 :    
32 :     template <typename REAL>
33 :     __details::image1d<REAL>::~image1d ()
34 :     {
35 :     free (this->_data);
36 :     }
37 :    
38 :     template <typename REAL>
39 :     bool image1d<REAL>::load (struct WorldBase *wrld, std::string const &name)
40 :     {
41 :     Nrrd *nin = load_nrrd_file (wrld, name);
42 :     if (nin == nullptr)
43 :     return true;
44 :     else {
45 :     bool sts = this->load (wrld, nin);
46 :     nrrdNix (nin);
47 :     return sts;
48 :     }
49 :     }
50 :    
51 :     template <typename REAL>
52 :     bool image1d<REAL>::load (struct WorldBase *wrld, const Nrrd *nin)
53 :     {
54 :     this->_img = std::make_shared<__details::image1d<REAL>>();
55 :    
56 :     // compute the offset to the first space axis
57 :     int base = nin->dim - nin->spaceDim;
58 :    
59 :     // FIXME: we should also be checking the shape of the voxels and the type of the samples.
60 :     // check that the image format is what we expect
61 :     if (nin->spaceDim != 1) {
62 :     error (wrld, "nrrd has unexpected dimension %d, expected 1\n", nin->spaceDim);
63 :     return true;
64 :     }
65 :     else if (nin->axis[base].kind != nrrdKindSpace) {
66 :     error (wrld, "nrrd has unexpected axis structure: %s\n",
67 :     airEnumStr(nrrdKind, nin->axis[base].kind));
68 :     return true;
69 :     }
70 :    
71 :     this->_img->_dim = 1;
72 :     this->_img->_size[0] = nin->axis[base+0].size;
73 :     this->_img->_data = nin->data;
74 :     this->_img->_dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin);
75 :    
76 :     //printf("SetImage \"%s\": space dim = 1, axes = <%d>\n",
77 :     //name, img->size[base+0]);
78 :    
79 :     // from the Nrrd file, we load the scaling and translation
80 :     double s, t;
81 :    
82 :     load_transform_1d (nin, s, t);
83 :    
84 :     this->_img->_s = static_cast<REAL>(1.0 / s);
85 :     this->_img->_t = static_cast<REAL>(-t / s);
86 :    
87 :     return false;
88 :     }
89 :    
90 :     /******************** 2D images ********************/
91 :    
92 :     static void load_transform_2d (const Nrrd *nin, mat2x2<double> &m, vec2<double> &t)
93 :     {
94 :     #define M(i,j) m[(i)*2+(j)]
95 :    
96 :     // compute the offset to the first space axis
97 :     int base = nin->dim - nin->spaceDim;
98 :    
99 :     // Image axis Scaling and Rotation
100 :     M(0,0) = nin->axis[base+0].spaceDirection[0];
101 :     M(1,0) = nin->axis[base+0].spaceDirection[1];
102 :     M(0,1) = nin->axis[base+1].spaceDirection[0];
103 :     M(1,1) = nin->axis[base+1].spaceDirection[1];
104 :    
105 :     // Image location
106 :     t[0] = nin->spaceOrigin[0];
107 :     t[1] = nin->spaceOrigin[1];
108 :    
109 :     #undef M
110 :     }
111 :    
112 :     inline double det2 (double a, double b, double c, double d)
113 :     {
114 :     return (a*d - b*c);
115 :     }
116 :    
117 :     /*! \brief compute the inverse of \arg m, storing the result in \arg i.
118 :     * \param m the matrix to invert
119 :     * \param i the inverted matrix
120 :     */
121 :     inline void invert (mat2x2<double> &i, mat2x2<double> &m)
122 :     {
123 :     double scale = 1.0 / det2(m[0], m[1], m[2], m[3]);
124 :    
125 :     i[0] = scale * m[3];
126 :     i[1] = -scale * m[1];
127 :     i[2] = -scale * m[2];
128 :     i[3] = scale * m[0];
129 :     }
130 :    
131 :     template <typename REAL>
132 :     __details::image2d<REAL>::~image2d ()
133 :     {
134 :     free (this->_data);
135 :     }
136 :    
137 :     template <typename REAL>
138 :     bool image2d<REAL>::load (struct WorldBase *wrld, std::string const &name)
139 :     {
140 :     Nrrd *nin = load_nrrd_file (wrld, name);
141 :     if (nin == nullptr)
142 :     return true;
143 :     else {
144 :     bool sts = this->load (wrld, nin);
145 :     nrrdNix (nin);
146 :     return sts;
147 :     }
148 :     }
149 :    
150 :     template <typename REAL>
151 :     bool image2d<REAL>::load (struct WorldBase *wrld, const Nrrd *nin)
152 :     {
153 :     this->_img = std::make_shared<__details::image2d<REAL>>();
154 :    
155 :     // compute the offset to the first space axis
156 :     int base = nin->dim - nin->spaceDim;
157 :    
158 :     // check that the image format is what we expect
159 :     if (nin->spaceDim != 2) {
160 :     error (wrld, "nrrd has unexpected dimension %d, expected 2\n", nin->spaceDim);
161 :     return true;
162 :     }
163 :     else if ((nin->axis[base].kind != nrrdKindSpace)
164 :     || (nin->axis[base+1].kind != nrrdKindSpace)) {
165 :     error (wrld, "nrrd has unexpected axis structure: %s %s\n",
166 :     airEnumStr(nrrdKind, nin->axis[base].kind),
167 :     airEnumStr(nrrdKind, nin->axis[base+1].kind));
168 :     return true;
169 :     }
170 :    
171 :     this->_img->_dim = 2;
172 :     this->_img->_size[0] = nin->axis[base+0].size;
173 :     this->_img->_size[1] = nin->axis[base+1].size;
174 :     this->_img->_data = nin->data;
175 :     this->_img->_dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin);
176 :    
177 :     // from the Nrrd file, we load the affine image-to-world transform matrix
178 :     mat2x2<double> m; // rotation and scaling
179 :     vec2<double> t; // translation part of the transform
180 :     load_transform_2d (nin, m, t);
181 :    
182 :     // compute inverse of m, which is the transform from the world basis to the image basis
183 :     mat2x2<double> mInv;
184 :     invert (mInv, m);
185 :    
186 :     // copy inverted matrix into the image data structure
187 :     for (int i = 0; i < 2; i++) {
188 :     for (int j = 0; j < 2; j++) {
189 :     REAL r = static_cast<REAL>(mInv[2*i+j]);
190 :     this->_img->_w2i[2*i+j] = r;
191 :     this->_img->_w2iT[2*j+i] = r; // transpose
192 :     }
193 :     }
194 :    
195 :     // transform the translation vector: inv([M t]) = [inv(M) -inv(M)t]
196 :     for (int i = 0; i < 2; i++) {
197 :     double sum = 0;
198 :     for (int j = 0; j < 2; j++) {
199 :     sum += mInv[2*i+j] * t[j];
200 :     }
201 :     this->_img->_tVec[i] = static_cast<REAL>(-sum);
202 :     }
203 :    
204 :     return false;
205 :     }
206 :    
207 :     /******************** 3D images ********************/
208 :    
209 :     static void load_transform_3d (const Nrrd *nin, mat3x3<double> &m, vec3<double> &t)
210 :     {
211 :     #define M(i,j) m[3*(i)+(j)]
212 :    
213 :     // compute the offset to the first space axis
214 :     int base = nin->dim - nin->spaceDim;
215 :    
216 :     // Image axis Scaling and Rotation
217 :     M(0,0) = nin->axis[base+0].spaceDirection[0];
218 :     M(1,0) = nin->axis[base+0].spaceDirection[1];
219 :     M(2,0) = nin->axis[base+0].spaceDirection[2];
220 :     M(0,1) = nin->axis[base+1].spaceDirection[0];
221 :     M(1,1) = nin->axis[base+1].spaceDirection[1];
222 :     M(2,1) = nin->axis[base+1].spaceDirection[2];
223 :     M(0,2) = nin->axis[base+2].spaceDirection[0];
224 :     M(1,2) = nin->axis[base+2].spaceDirection[1];
225 :     M(2,2) = nin->axis[base+2].spaceDirection[2];
226 :    
227 :     // Image location
228 :     t[0] = nin->spaceOrigin[0];
229 :     t[1] = nin->spaceOrigin[1];
230 :     t[2] = nin->spaceOrigin[2];
231 :    
232 :     #undef M
233 :     }
234 :    
235 :     /*! \brief compute the inverse of \arg m, storing the result in \arg i.
236 :     * \param m the matrix to invert
237 :     * \param i the inverted matrix
238 :     */
239 :     inline void invert (mat3x3<double> &i, mat3x3<double> &m)
240 :     {
241 :     double det =
242 :     ( m[0]*m[4]*m[8]
243 :     + m[3]*m[7]*m[2]
244 :     + m[6]*m[1]*m[5]
245 :     - m[6]*m[4]*m[2]
246 :     - m[3]*m[1]*m[8]
247 :     - m[0]*m[7]*m[5]);
248 :     double scale = 1.0 / det;
249 :    
250 :     i[0] = scale * det2 (m[4], m[5], m[7], m[8]);
251 :     i[1] = scale * det2 (m[2], m[1], m[8], m[7]);
252 :     i[2] = scale * det2 (m[1], m[2], m[4], m[5]);
253 :     i[3] = scale * det2 (m[5], m[3], m[8], m[6]);
254 :     i[4] = scale * det2 (m[0], m[2], m[6], m[8]);
255 :     i[5] = scale * det2 (m[2], m[0], m[5], m[3]);
256 :     i[6] = scale * det2 (m[3], m[4], m[6], m[7]);
257 :     i[7] = scale * det2 (m[1], m[0], m[7], m[6]);
258 :     i[8] = scale * det2 (m[0], m[1], m[3], m[4]);
259 :     }
260 :    
261 :     template <typename REAL>
262 :     __details::image3d<REAL>::~image3d ()
263 :     {
264 :     free (this->_data);
265 :     }
266 :    
267 :     template <typename REAL>
268 :     bool image3d<REAL>::load (struct WorldBase *wrld, std::string const &name)
269 :     {
270 :     Nrrd *nin = load_nrrd_file (wrld, name);
271 :     if (nin == nullptr)
272 :     return true;
273 :     else {
274 :     bool sts = this->load (wrld, nin);
275 :     nrrdNix (nin);
276 :     return sts;
277 :     }
278 :     }
279 :    
280 :     template <typename REAL>
281 :     bool image3d<REAL>::load (struct WorldBase *wrld, const Nrrd *nin)
282 :     {
283 :     this->_img = std::make_shared<__details::image3d<REAL>>();
284 :    
285 :     // compute the offset to the first space axis
286 :     int base = nin->dim - nin->spaceDim;
287 :    
288 :     // check that the image format is what we expect
289 :     if (nin->spaceDim != 3) {
290 :     error (wrld, "nrrd has unexpected dimension %d, expected 3\n", nin->spaceDim);
291 :     return true;
292 :     }
293 :     else if ((nin->axis[base].kind != nrrdKindSpace)
294 :     || (nin->axis[base+1].kind != nrrdKindSpace)
295 :     || (nin->axis[base+2].kind != nrrdKindSpace)) {
296 :     error (wrld, "nrrd has unexpected axis structure: %s %s %s\n",
297 :     airEnumStr(nrrdKind, nin->axis[base].kind),
298 :     airEnumStr(nrrdKind, nin->axis[base+1].kind),
299 :     airEnumStr(nrrdKind, nin->axis[base+2].kind));
300 :     return true;
301 :     }
302 :    
303 :     this->_img->_dim = 3;
304 :     this->_img->_size[0] = nin->axis[base+0].size;
305 :     this->_img->_size[1] = nin->axis[base+1].size;
306 :     this->_img->_size[2] = nin->axis[base+2].size;
307 :     this->_img->_data = nin->data;
308 :     this->_img->_dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin);
309 :    
310 :     // from the Nrrd file, we load the affine image-to-world transform matrix
311 :     mat3x3<double> m; // rotation and scaling
312 :     vec3<double> t; // translation part of the transform
313 :     load_transform_3d (nin, m, t);
314 :    
315 :     // compute inverse of m, which is the transform from the world basis to the image basis
316 :     mat3x3<double> mInv;
317 :     invert (mInv, m);
318 :    
319 :     // copy results into the image data structure
320 :     for (int i = 0; i < 3; i++) {
321 :     for (int j = 0; j < 3; j++) {
322 :     REAL r = static_cast<REAL>(mInv[2*i+j]);
323 :     this->_img->_w2i[3*i+j] = r;
324 :     this->_img->_w2iT[3*j+i] = r; // transpose
325 :     }
326 :     }
327 :    
328 :     // transform the translation vector: inv([M t]) = [inv(M) -inv(M)t]
329 :     for (int i = 0; i < 3; i++) {
330 :     double sum = 0;
331 :     for (int j = 0; j < 3; j++) {
332 :     sum += mInv[3*i+j] * t[j];
333 :     }
334 :     this->_img->_tVec[i] = static_cast<REAL>(-sum);
335 :     }
336 :    
337 :     return false;
338 :     }
339 :    
340 :     } // namespace Diderot
341 :    
342 :     void foo (std::string name)
343 :     {
344 :     Diderot::image1d<float> img1;
345 :     Diderot::image2d<float> img2;
346 :     Diderot::image3d<float> img3;
347 :    
348 :     img1.load(nullptr, name);
349 :     img2.load(nullptr, name);
350 :     img3.load(nullptr, name);
351 :     }

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