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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : jhr 5597 /*! \file image.cxx
2 :     *
3 :     * Utility code to support images in Diderot. Most of the actual implementation can be
4 :     * found in diderot/image-inst.hxx; this file just contains the shared non-template code.
5 :     *
6 :     * \author John Reppy
7 :     */
8 :    
9 :     /*
10 :     * This code is part of the Diderot Project (http://diderot-language.cs.uchicago.edu)
11 :     *
12 :     * COPYRIGHT (c) 2016 The University of Chicago
13 :     * All rights reserved.
14 :     */
15 :    
16 :     #include <cstdlib>
17 : adrianlehm 5619 #include "diderot/image.h"
18 : jhr 5597
19 :     namespace diderot {
20 :    
21 :     namespace __details {
22 :    
23 :     void load_transform_1d (const Nrrd *nin, double &s, double &t)
24 :     {
25 :     // compute the offset to the first space axis
26 :     int base = nin->dim - nin->spaceDim;
27 :    
28 :     // Image axis Scaling and Rotation
29 :     s = nin->axis[base+0].spaceDirection[0];
30 :    
31 :     // Image location
32 :     t = nin->spaceOrigin[0];
33 :     }
34 :    
35 :    
36 :     void load_transform_2d (const Nrrd *nin, mat2x2<double> &m, real2<double> &t)
37 :     {
38 :     #define M(i,j) m[(i)*2+(j)]
39 :    
40 :     // compute the offset to the first space axis
41 :     int base = nin->dim - nin->spaceDim;
42 :    
43 :     // Image axis Scaling and Rotation
44 :     M(0,0) = nin->axis[base+0].spaceDirection[0];
45 :     M(1,0) = nin->axis[base+0].spaceDirection[1];
46 :     M(0,1) = nin->axis[base+1].spaceDirection[0];
47 :     M(1,1) = nin->axis[base+1].spaceDirection[1];
48 :    
49 :     // Image location
50 :     t[0] = nin->spaceOrigin[0];
51 :     t[1] = nin->spaceOrigin[1];
52 :    
53 :     #undef M
54 :     }
55 :    
56 :     inline double det2 (double a, double b, double c, double d)
57 :     {
58 :     return (a*d - b*c);
59 :     }
60 :    
61 :     void invert2x2 (mat2x2<double> &i, mat2x2<double> &m)
62 :     {
63 :     double scale = 1.0 / det2(m[0], m[1], m[2], m[3]);
64 :    
65 :     i[0] = scale * m[3];
66 :     i[1] = -scale * m[1];
67 :     i[2] = -scale * m[2];
68 :     i[3] = scale * m[0];
69 :     }
70 :    
71 :     void load_transform_3d (const Nrrd *nin, mat3x3<double> &m, real3<double> &t)
72 :     {
73 :     #define M(i,j) m[3*(i)+(j)]
74 :    
75 :     // compute the offset to the first space axis
76 :     int base = nin->dim - nin->spaceDim;
77 :    
78 :     // Image axis Scaling and Rotation
79 :     M(0,0) = nin->axis[base+0].spaceDirection[0];
80 :     M(1,0) = nin->axis[base+0].spaceDirection[1];
81 :     M(2,0) = nin->axis[base+0].spaceDirection[2];
82 :     M(0,1) = nin->axis[base+1].spaceDirection[0];
83 :     M(1,1) = nin->axis[base+1].spaceDirection[1];
84 :     M(2,1) = nin->axis[base+1].spaceDirection[2];
85 :     M(0,2) = nin->axis[base+2].spaceDirection[0];
86 :     M(1,2) = nin->axis[base+2].spaceDirection[1];
87 :     M(2,2) = nin->axis[base+2].spaceDirection[2];
88 :    
89 :     // Image location
90 :     t[0] = nin->spaceOrigin[0];
91 :     t[1] = nin->spaceOrigin[1];
92 :     t[2] = nin->spaceOrigin[2];
93 :    
94 :     #undef M
95 :     }
96 :    
97 :     /*! \brief compute the inverse of \arg m, storing the result in \arg i.
98 :     * \param m the matrix to invert
99 :     * \param i the inverted matrix
100 :     */
101 :     void invert3x3 (mat3x3<double> &i, mat3x3<double> &m)
102 :     {
103 :     double det =
104 :     ( m[0]*m[4]*m[8]
105 :     + m[3]*m[7]*m[2]
106 :     + m[6]*m[1]*m[5]
107 :     - m[6]*m[4]*m[2]
108 :     - m[3]*m[1]*m[8]
109 :     - m[0]*m[7]*m[5]);
110 :     double scale = 1.0 / det;
111 :    
112 :     i[0] = scale * det2 (m[4], m[5], m[7], m[8]);
113 :     i[1] = scale * det2 (m[2], m[1], m[8], m[7]);
114 :     i[2] = scale * det2 (m[1], m[2], m[4], m[5]);
115 :     i[3] = scale * det2 (m[5], m[3], m[8], m[6]);
116 :     i[4] = scale * det2 (m[0], m[2], m[6], m[8]);
117 :     i[5] = scale * det2 (m[2], m[0], m[5], m[3]);
118 :     i[6] = scale * det2 (m[3], m[4], m[6], m[7]);
119 :     i[7] = scale * det2 (m[1], m[0], m[7], m[6]);
120 :     i[8] = scale * det2 (m[0], m[1], m[3], m[4]);
121 :     }
122 :    
123 :     } // namespace __details
124 :    
125 :     /******************** nrrd_proxy ********************/
126 :    
127 :     bool nrrd_proxy::check_nrrd (struct world_base *wrld, const Nrrd *nin)
128 :     {
129 :     // compute the offset to the first space axis
130 :     int base = nin->dim - nin->spaceDim;
131 :    
132 :     // check that the image format is what we expect
133 :     if (nin->spaceDim != this->_dim) {
134 :     wrld->error ("nrrd has unexpected dimension %d, expected %d\n",
135 :     nin->spaceDim, this->_dim);
136 :     return true;
137 :     }
138 :    
139 :     // check the number of elements per voxel
140 :     int nElems = 1;
141 :     for (int i = 0; i < base; i++) {
142 :     nElems *= nin->axis[i].size;
143 :     }
144 :     if (nElems != this->_voxelSz) {
145 :     wrld->error ("nrrd has %d elements per voxel (expected %d)\n",
146 :     nElems, this->_voxelSz);
147 :     return true;
148 :     }
149 :    
150 :     // check the axis kinds and sizes
151 :     bool bad = false;
152 :     for (int i = 0; i < this->_dim; i++) {
153 :     if (nin->axis[base+i].kind != nrrdKindSpace) {
154 :     wrld->error ("nrrd axis[%d] has kind %s (expected %s)\n",
155 :     i, airEnumStr(nrrdKind, nin->axis[base+i].kind),
156 :     airEnumStr(nrrdKind, nrrdKindSpace));
157 :     bad = true;
158 :     }
159 :     if (nin->axis[base+i].size != this->_sizes[i]) {
160 :     wrld->error ("nrrd axis[%d] has size %d (expected %d)\n",
161 :     i, nin->axis[base+i].size, this->_sizes[i]);
162 :     bad = true;
163 :     }
164 :     }
165 :    
166 :     return bad;
167 :     }
168 :    
169 :     } // namespace Diderot

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