EMAN2
marchingcubes.cpp
Go to the documentation of this file.
1/*
2 * Author: Tao Ju, 5/16/2007 <taoju@cs.wustl.edu>, code ported by Grant Tang
3 * code extensively modified and optimized by David Woolford
4 * Copyright (c) 2000-2006 Baylor College of Medicine
5 *
6 * This software is issued under a joint BSD/GNU license. You may use the
7 * source code in this file under either license. However, note that the
8 * complete EMAN2 and SPARX software packages have some GPL dependencies,
9 * so you are responsible for compliance with the licenses of these packages
10 * if you opt to use BSD licensing. The warranty disclaimer below holds
11 * in either instance.
12 *
13 * This complete copyright notice must be included in any revised version of the
14 * source code. Additional authorship citations may be added, but existing
15 * author citations must be preserved.
16 *
17 * This program is free software; you can redistribute it and/or modify
18 * it under the terms of the GNU General Public License as published by
19 * the Free Software Foundation; either version 2 of the License, or
20 * (at your option) any later version.
21 *
22 * This program is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with this program; if not, write to the Free Software
29 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30 *
31 * */
32#ifdef USE_OPENGL
33
34#ifdef _WIN32
35 #include <windows.h>
36#endif //_WIN32
37
38#ifndef GL_GLEXT_PROTOTYPES
39 #define GL_GLEXT_PROTOTYPES
40#endif //GL_GLEXT_PROTOTYPES
41
42
43#include "marchingcubes.h"
44
45#include <time.h>
46#include <math.h>
47
48using namespace EMAN;
49#include "transform.h"
50#include "emobject.h"
51#include "vec3.h"
52
53//a2fVertexOffset lists the positions, relative to vertex0, of each of the 8 vertices of a cube
54static const int a2fVertexOffset[8][3] =
55{
56 {0, 0, 0},{1, 0, 0},{1, 1, 0},{0, 1, 0},
57 {0, 0, 1},{1, 0, 1},{1, 1, 1},{0, 1, 1}
58};
59
60
61static const int a2fPosXOffset[4][3] =
62{
63 {2, 0, 0},{2, 1, 0},
64 {2, 0, 1},{2, 1, 1}
65};
66
67static const int a2fPosYOffset[4][3] =
68{
69 {1, 2, 0},{0, 2, 0},
70 {1, 2, 1},{0, 2, 1}
71};
72
73static const int a2fPosZOffset[4][3] =
74{
75 {0, 0, 2},{1, 0, 2},
76 {1, 1, 2},{0, 1, 2}
77};
78
79static const int a2fPosXPosYOffset[2][3] =
80{
81 {2, 2, 0},{2, 2, 1}
82};
83
84static const int a2fPosXPosZOffset[2][3] =
85{
86 {2, 0, 2},{2, 1, 2}
87};
88
89static const int a2fPosYPosZOffset[2][3] =
90{
91 {1, 2, 2},{0, 2, 2}
92};
93
94
95static const int a2fPosXPosZPosYOffset[3] =
96{
97 2, 2, 2
98};
99
100//a2fVertexOffset lists the positions, relative to vertex0, of each of the 8 vertices of a cube
101static const int a2OddXOffset[8] =
102{
103 1,0,0,1,
104 1,0,0,1
105};
106
107static const int a2OddYOffset[8] =
108{
109 1,1,0,0,
110 1,1,0,0
111};
112
113static const int a2OddZOffset[8] =
114{
115 1,1,1,1,
116 0,0,0,0
117};
118
119//a2iEdgeConnection lists the index of the endpoint vertices for each of the 12 edges of the cube
120static const int a2iEdgeConnection[12][2] =
121{
122 {0,1}, {1,2}, {2,3}, {3,0},
123 {4,5}, {5,6}, {6,7}, {7,4},
124 {0,4}, {1,5}, {2,6}, {3,7}
125};
126
127//a2fEdgeDirection lists the direction vector (vertex1-vertex0) for each edge in the cube
128static const float a2fEdgeDirection[12][3] =
129{
130 {1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{-1.0, 0.0, 0.0},{0.0, -1.0, 0.0},
131 {1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{-1.0, 0.0, 0.0},{0.0, -1.0, 0.0},
132 {0.0, 0.0, 1.0},{0.0, 0.0, 1.0},{ 0.0, 0.0, 1.0},{0.0, 0.0, 1.0}
133};
134
135// right = 0, up = 1, back = 2
136static const int edgeLookUp[12][4] =
137{
138 {0,0,0,0},{1,0,0,1},{0,1,0,0},{0,0,0,1},
139 {0,0,1,0},{1,0,1,1},{0,1,1,0},{0,0,1,1},
140 {0,0,0,2},{1,0,0,2},{1,1,0,2},{0,1,0,2}
141};
142
143// For any edge, if one vertex is inside of the surface and the other is outside of the surface
144// then the edge intersects the surface
145// For each of the 8 vertices of the cube can be two possible states : either inside or outside of the surface
146// For any cube the are 2^8=256 possible sets of vertex states
147// This table lists the edges intersected by the surface for all 256 possible vertex states
148// There are 12 edges. For each entry in the table, if edge #n is intersected, then bit #n is set to 1
149
150int aiCubeEdgeFlags[256]=
151{
152 0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
153 0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
154 0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435, 0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
155 0x3a0, 0x2a9, 0x1a3, 0x0aa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
156 0x460, 0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
157 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
158 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c, 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
159 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0x0cc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
160 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0x0cc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
161 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x055, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
162 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5, 0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
163 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460,
164 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0,
165 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x033, 0x339, 0x230,
166 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190,
167 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x000
168};
169
170// For each of the possible vertex states listed in aiCubeEdgeFlags there is a specific triangulation
171// of the edge intersection points. a2iTriangleConnectionTable lists all of thminvals[cur_level]-> in the form of
172// 0-5 edge triples with the list terminated by the invalid value -1.
173// For example: a2iTriangleConnectionTable[3] list the 2 triangles formed when corner[0]
174// and corner[1] are inside of the surface, but the rest of the cube is not.
175//
176// I found this table in an example program someone wrote long ago. It was probably generated by hand
177
178int a2iTriangleConnectionTable[256][16] =
179{
180 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
181 {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
182 {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
183 {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
184 {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
185 {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
186 {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
187 {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
188 {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
189 {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
190 {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
191 {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
192 {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
193 {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
194 {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
195 {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
196 {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
197 {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
198 {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
199 {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
200 {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
201 {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
202 {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
203 {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
204 {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
205 {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
206 {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
207 {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
208 {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
209 {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
210 {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
211 {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
212 {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
213 {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
214 {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
215 {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
216 {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
217 {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
218 {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
219 {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
220 {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
221 {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
222 {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
223 {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
224 {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
225 {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
226 {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
227 {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
228 {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
229 {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
230 {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
231 {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
232 {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
233 {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
234 {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
235 {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
236 {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
237 {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
238 {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
239 {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
240 {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
241 {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
242 {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
243 {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
244 {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
245 {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
246 {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
247 {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
248 {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
249 {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
250 {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
251 {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
252 {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
253 {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
254 {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
255 {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
256 {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
257 {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
258 {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
259 {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
260 {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
261 {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
262 {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
263 {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
264 {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
265 {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
266 {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
267 {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
268 {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
269 {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
270 {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
271 {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
272 {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
273 {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
274 {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
275 {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
276 {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
277 {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
278 {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
279 {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
280 {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
281 {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
282 {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
283 {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
284 {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
285 {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
286 {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
287 {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
288 {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
289 {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
290 {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
291 {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
292 {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
293 {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
294 {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
295 {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
296 {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
297 {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
298 {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
299 {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
300 {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
301 {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
302 {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
303 {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
304 {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
305 {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
306 {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
307 {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
308 {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
309 {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
310 {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
311 {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
312 {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
313 {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
314 {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
315 {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
316 {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
317 {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
318 {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
319 {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
320 {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
321 {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
322 {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
323 {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
324 {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
325 {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
326 {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
327 {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
328 {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
329 {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
330 {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
331 {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
332 {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
333 {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
334 {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
335 {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
336 {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
337 {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
338 {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
339 {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
340 {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
341 {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
342 {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
343 {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
344 {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
345 {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
346 {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
347 {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
348 {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
349 {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
350 {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
351 {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
352 {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
353 {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
354 {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
355 {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
356 {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
357 {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
358 {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
359 {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
360 {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
361 {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
362 {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
363 {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
364 {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
365 {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
366 {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
367 {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
368 {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
369 {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
370 {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
371 {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
372 {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
373 {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
374 {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
375 {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
376 {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
377 {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
378 {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
379 {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
380 {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
381 {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
382 {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
383 {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
384 {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
385 {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
386 {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
387 {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
388 {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
389 {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
390 {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
391 {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
392 {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
393 {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
394 {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
395 {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
396 {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
397 {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
398 {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
399 {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
400 {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
401 {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
402 {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
403 {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
404 {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
405 {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
406 {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
407 {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
408 {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
409 {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
410 {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
411 {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
412 {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
413 {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
414 {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
415 {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
416 {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
417 {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
418 {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
419 {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
420 {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
421 {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
422 {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
423 {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
424 {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
425 {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
426 {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
427 {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
428 {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
429 {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
430 {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
431 {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
432 {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
433 {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
434 {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
435 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
436};
437
438ColorRGBGenerator::ColorRGBGenerator()
439 : rgbmode(0), originx(0), originy(0), originz(0), inner(0.0), outer(0.0), minimum(0.0), maximum(0.0), needtorecolor(1), colormap(0), em_data(0)
440{
441}
442
443ColorRGBGenerator::ColorRGBGenerator(EMData* data)
444 : rgbmode(0), minimum(0.0), maximum(0.0), needtorecolor(1), colormap(0)
445{
446 set_data(data);
447}
448
449void ColorRGBGenerator::set_data(EMData* data)
450{
451 em_data = data;
452 originx = data->get_xsize()/2;
453 originy = data->get_ysize()/2;
454 originz = data->get_zsize()/2;
455 inner = 0;
456 outer = (float)originx;
457}
458
459void ColorRGBGenerator::set_cmap_data(EMData* data)
460{
461 cmap = data;
462 setMinMax(data->get_attr("minimum"), data->get_attr("maximum"));
463}
464
465void ColorRGBGenerator::generateRadialColorMap()
466{
467
468 int size = int(sqrt((float)originx*originx + originy*originy + originz*originz));
469
470 if(colormap) delete colormap;
471 colormap = new float[size*3];
472
473 for(int i = 0; i < size; i++){
474 float normrad = 4.189f*(i - inner)/(outer - inner);
475 if(normrad < 2.094){
476 if (normrad < 0.0) normrad = 0.0;
477 colormap[i*3] = 0.5f*(1 + cos(normrad)/cos(1.047f - normrad));
478 colormap[i*3 + 1] = 1.5f - colormap[i*3];
479 colormap[i*3 + 2] = 0.0;
480 }
481 if(normrad >= 2.094){
482 if (normrad > 4.189f) normrad = 4.189f;
483 normrad -= 2.094f;
484 colormap[i*3] = 0.0;
485 colormap[i*3 + 1] = 0.5f*(1 + cos(normrad)/cos(1.047f - normrad));
486 colormap[i*3 + 2] = 1.5f - colormap[i*3 + 1];
487 }
488 }
489}
490
491float* ColorRGBGenerator::getRGBColor(int x, int y, int z)
492{
493 const float brt=0.7;
494 // Get data using radial color info
495 if (rgbmode == 1){
496 //calculate radius
497 float rad = sqrt(float(pow(x-originx,2) + pow(y-originy,2) + pow(z-originz,2)));
498
499 //This indicates that new color info needs to be bound (to the GPU)
500 if(needtorecolor){
501 generateRadialColorMap();
502 needtorecolor = false;
503 }
504
505 return &colormap[int(rad)*3];
506 }
507 // Get data using color map
508 if (rgbmode == 2){
509 float value = cmap->get_value_at(x, y, z);
510 value = 4.189f*(value - minimum)/(maximum - minimum);
511 if(value < 2.094){
512 if (value < 0.0) value = 0.0;
513 rgb[0] = brt* 0.5f*(1 + cos(value)/cos(1.047f - value));
514 rgb[1] = brt* (1.5f - rgb[0]);
515 rgb[2] = 0.0;
516 }
517 if(value >= 2.094){
518 if (value > 4.189f) value = 4.189f;
519 value -= 2.094f;
520 rgb[0] = 0.0;
521 rgb[1] = brt * 0.5f*(1 + cos(value)/cos(1.047f - value));
522 rgb[2] = brt * (1.5f - rgb[1]);
523 }
524
525 return &rgb[0];
526 }
527
528 return &colormap[0];
529}
530
531MarchingCubes::MarchingCubes()
532 : _isodl(0), needtobind(1)
533{
534
535if ((int(glGetString(GL_VERSION)[0])-48)>2){
536 rgbgenerator = ColorRGBGenerator();
537
538// #ifdef _WIN32
539// typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
540// PFNGLGENBUFFERSPROC glGenBuffers;
541// glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
542// #endif //_WIN32
543//
544// glGenBuffers(4, buffer);
545}
546
547}
548
549MarchingCubes::MarchingCubes(EMData * em)
550 : _isodl(0)
551{
552if ((int(glGetString(GL_VERSION)[0])-48)>2){
553 rgbgenerator = ColorRGBGenerator();
554
555// #ifdef _WIN32
556// typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
557// PFNGLGENBUFFERSPROC glGenBuffers;
558// glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
559// #endif //_WIN32
560//
561// glGenBuffers(4, buffer);
562 set_data(em);}
563else{
564 set_data(em);
565 }
566}
567
568
569
570void MarchingCubes::clear_min_max_vals()
571{
572 for (vector<EMData*>::iterator it = minvals.begin(); it != minvals.end(); ++it)
573 {
574 if ( (*it) != 0 ) delete *it;
575 }
576 minvals.clear();
577
578 for (vector<EMData*>::iterator it = maxvals.begin(); it != maxvals.end(); ++it)
579 {
580 if ( (*it) != 0 ) delete *it;
581 }
582 maxvals.clear();
583}
584
585bool MarchingCubes::calculate_min_max_vals()
586{
587
588 if (_emdata == NULL ) throw NullPointerException("Error, cannot generate search tree if the overriding EMData object is NULL");
589
590 clear_min_max_vals();
591
592 int nx = _emdata->get_xsize();
593 int ny = _emdata->get_ysize();
594 int nz = _emdata->get_zsize();
595
596 // Create the binary min max tree
597 while ( nx > 1 || ny > 1 || nz > 1 )
598 {
599 int size = minvals.size();
600
601 if ( size == 0 ){
602 Dict a;
603 // Setting search to 3 at the bottom level is most important.
604 // FIXME - d.woolford add comments heere
605 a["search"] = 3;
606 minvals.push_back(_emdata->process("math.minshrink",a));
607 maxvals.push_back(_emdata->process("math.maxshrink",a));
608 }else {
609 minvals.push_back(minvals[size-1]->process("math.minshrink"));
610 maxvals.push_back(maxvals[size-1]->process("math.maxshrink"));
611 }
612
613 nx = minvals[size]->get_xsize();
614 ny = minvals[size]->get_ysize();
615 nz = minvals[size]->get_zsize();
616#if MARCHING_CUBES_DEBUG
617 cout << "dims are " << nx << " " << ny << " " << nz << endl;
618#endif
619 }
620
621 drawing_level = -1;
622
623 return true;
624}
625
626MarchingCubes::~MarchingCubes() {
627 clear_min_max_vals();
628
629if ((int(glGetString(GL_VERSION)[0])-48)>2){
630// #ifdef _WIN32
631// typedef void (APIENTRYP PFNGLDELETEBUFFERSPROC) (GLsizei n, const GLuint *buffers);
632// PFNGLDELETEBUFFERSPROC glDeleteBuffers;
633// glDeleteBuffers = (PFNGLDELETEBUFFERSPROC) wglGetProcAddress("glDeleteBuffers");
634// #endif //_WIN32
635//
636// glDeleteBuffers(4, buffer);
637}
638}
639
640Dict MarchingCubes::get_isosurface()
641{
642 calculate_surface();
643 Dict d;
644 d.put("points", (float*)pp.get_data());
645 for (unsigned int i = 0; i < ff.elem(); ++i ) ff[i] /= 3;
646 d.put("faces", (unsigned int*)ff.get_data());
647 d.put("normals", (float*)nn.get_data());
648 d.put("size", ff.elem());
649 return d;
650}
651
652void MarchingCubes::surface_face_z()
653{
654 float* f = pp.get_data();
655 float* n = nn.get_data();
656 for (unsigned int i = 0; i < pp.elem(); i += 3 ) {
657 if (f[i+2] == 0.5) continue;
658 Vec3f z(0,0,1.0);
659 Vec3f axis(-n[i+1],n[i],0);
660 axis.normalize();
661
662 Dict d;
663 d["type"] = "spin";
664 d["Omega"] = 90.f;
665 d["n1"] = axis[0];
666 d["n2"] = axis[1];
667 d["n3"] = 0;
668 Transform t(d);
669 Vec3f delta = t*z;
670
671 f[i] += delta[0]*.25f;
672 f[i+1] += delta[1]*.25f;
673 f[i+2] = 0.5;
674 }
675
676 for (unsigned int i = 0; i < nn.elem(); i += 3 ) {
677 n[i] = 0;
678 n[i+1] = 0;
679 n[i+2] = 1;
680 }
681}
682
683void MarchingCubes::set_data(EMData* data)
684{
685 if ( data->get_zsize() == 1 ) throw ImageDimensionException("The z dimension of the image must be greater than 1");
686 _emdata = data;
687 calculate_min_max_vals();
688 rgbgenerator.set_data(data);
689}
690
691void MarchingCubes::set_surface_value(const float value) {
692
693 if(_surf_value == value) return;
694
695 _surf_value = value;
696
697}
698
699void MarchingCubes::calculate_surface() {
700
701 if ( _emdata == 0 ) throw NullPointerException("Error, attempt to generate isosurface, but the emdata image object has not been set");
702 if ( minvals.size() == 0 || maxvals.size() == 0 ) throw NotExistingObjectException("Vector of EMData pointers", "Error, the min and max val search trees have not been created");
703
704 point_map.clear();
705 pp.clear();
706 nn.clear();
707 ff.clear();
708 vv.clear();
709
710#if MARCHING_CUBES_DEBUG
711 int time0 = clock();
712#endif
713
714 float min = minvals[minvals.size()-1]->get_value_at(0,0,0);
715 float max = maxvals[minvals.size()-1]->get_value_at(0,0,0);
716 if ( min < _surf_value && max > _surf_value) draw_cube(0,0,0,minvals.size()-1);
717
718#if MARCHING_CUBES_DEBUG
719 int time1 = clock();
720 cout << "It took " << (time1-time0) << " " << (float)(time1-time0)/CLOCKS_PER_SEC << " to traverse the search tree and generate polygons" << endl;
721 cout << "... using surface value " << _surf_value << endl;
722#endif
723}
724
725void MarchingCubes::draw_cube(const int x, const int y, const int z, const int cur_level ) {
726
727 if ( cur_level == drawing_level )
728 {
729 EMData* e;
730 if ( drawing_level == - 1 ) e = _emdata;
731 else e = minvals[drawing_level];
732 if ( x < (e->get_xsize()-1) && y < (e->get_ysize()-1) && z < (e->get_zsize()-1))
733 marching_cube(x,y,z, cur_level);
734 }
735 else
736 {
737 EMData* e;
738 if ( cur_level > 0 ) {
739 int xsize = minvals[cur_level-1]->get_xsize();
740 int ysize = minvals[cur_level-1]->get_ysize();
741 int zsize = minvals[cur_level-1]->get_zsize();
742 e = minvals[cur_level-1];
743 for(int i=0; i<8; ++i ) {
744 int xx = 2*x+a2fVertexOffset[i][0];
745 if ( xx >= xsize ) continue;
746 int yy = 2*y+a2fVertexOffset[i][1];
747 if ( yy >= ysize ) continue;
748 int zz = 2*z+a2fVertexOffset[i][2];
749 if ( zz >= zsize ) continue;
750
751 float min = minvals[cur_level-1]->get_value_at(xx,yy,zz);
752 float max = maxvals[cur_level-1]->get_value_at(xx,yy,zz);
753 if ( min < _surf_value && max > _surf_value)
754 draw_cube(xx,yy,zz,cur_level-1);
755 }
756 }
757 else {
758 e = _emdata;
759 for(int i=0; i<8; ++i ) {
760 draw_cube(2*x+a2fVertexOffset[i][0],2*y+a2fVertexOffset[i][1],2*z+a2fVertexOffset[i][2],cur_level-1);
761 }
762 }
763
764 if ( x == (minvals[cur_level]->get_xsize()-1) ) {
765 if ( e->get_xsize() > 2*x ){
766 for(int i=0; i<4; ++i ) {
767 draw_cube(2*x+a2fPosXOffset[i][0],2*y+a2fPosXOffset[i][1],2*z+a2fPosXOffset[i][2],cur_level-1);
768 }
769 }
770 if ( y == (minvals[cur_level]->get_ysize()-1) ) {
771 if ( e->get_ysize() > 2*y ) {
772 for(int i=0; i<2; ++i ) {
773 draw_cube(2*x+a2fPosXPosYOffset[i][0],2*y+a2fPosXPosYOffset[i][1],2*z+a2fPosXPosYOffset[i][2],cur_level-1);
774 }
775 }
776 if ( z == (minvals[cur_level]->get_zsize()-1) ){
777 if ( e->get_zsize() > 2*z ) {
778 draw_cube(2*x+2,2*y+2,2*z+2,cur_level-1);
779 }
780 }
781 }
782 if ( z == (minvals[cur_level]->get_zsize()-1) ) {
783 if ( e->get_zsize() > 2*z ) {
784 for(int i=0; i<2; ++i ) {
785 draw_cube(2*x+a2fPosXPosZOffset[i][0],2*y+a2fPosXPosZOffset[i][1],2*z+a2fPosXPosZOffset[i][2],cur_level-1);
786 }
787 }
788 }
789 }
790 if ( y == (minvals[cur_level]->get_ysize()-1) ) {
791 if ( e->get_ysize() > 2*y ) {
792 for(int i=0; i<4; ++i ) {
793 draw_cube(2*x+a2fPosYOffset[i][0],2*y+a2fPosYOffset[i][1],2*z+a2fPosYOffset[i][2],cur_level-1);
794 }
795 }
796 if ( z == (minvals[cur_level]->get_ysize()-1) ) {
797 if ( e->get_zsize() > 2*z ) {
798 for(int i=0; i<2; ++i ) {
799 draw_cube(2*x+a2fPosYPosZOffset[i][0],2*y+a2fPosYPosZOffset[i][1],2*z+a2fPosYPosZOffset[i][2],cur_level-1);
800 }
801 }
802 }
803 }
804 if ( z == (minvals[cur_level]->get_zsize()-1) ) {
805 if ( e->get_zsize() > 2*z ) {
806 for(int i=0; i<4; ++i ) {
807 draw_cube(2*x+a2fPosZOffset[i][0],2*y+a2fPosZOffset[i][1],2*z+a2fPosZOffset[i][2],cur_level-1);
808 }
809 }
810 }
811
812 }
813}
814
815void MarchingCubes::get_normal(Vector3 &normal, int fX, int fY, int fZ)
816{
817 normal[0] = _emdata->get_value_at(fX-1, fY, fZ) - _emdata->get_value_at(fX+1, fY, fZ);
818 normal[1] = _emdata->get_value_at(fX, fY-1, fZ) - _emdata->get_value_at(fX, fY+1, fZ);
819 normal[2] = _emdata->get_value_at(fX, fY, fZ-1) - _emdata->get_value_at(fX, fY, fZ+1);
820 normal.normalize();
821}
822
823float MarchingCubes::get_offset(float fValue1, float fValue2, float fValueDesired)
824{
825 float fDelta = fValue2 - fValue1;
826
827 if(fDelta == 0.0f)
828 {
829 return 0.5f;
830 }
831 return (fValueDesired - fValue1)/fDelta;
832}
833
834int MarchingCubes::get_edge_num(int x, int y, int z, int edge) {
835 // edge direction is right, down, back (x, y, z)
836 unsigned int index = 0;
837 index = (x << 22) | (y << 12) | (z << 2) | edge;
838 return index;
839}
840
841void MarchingCubes::color_vertices()
842{
843 cc.clear();
844 int scaling = pow(2,drawing_level + 1); // Needed to account for sampling rate
845 //Color vertices. We don't need to rerun marching cubes on color vertices, so this method improves effciency
846 for(unsigned int i = 0; i < vv.elem(); i+=3){
847 float* color = rgbgenerator.getRGBColor(scaling*vv[i], scaling*vv[i+1], scaling*vv[i+2]);
848 cc.push_back_3(color);
849 }
850 rgbgenerator.setNeedToRecolor(false);
851}
852
853void MarchingCubes::marching_cube(int fX, int fY, int fZ, int cur_level)
854{
855// extern int aiCubeEdgeFlags[256];
856// extern int a2iTriangleConnectionTable[256][16];
857
858 int iCorner, iVertex, iVertexTest, iEdge, iTriangle, iFlagIndex, iEdgeFlags;
859 float fOffset;
860 Vector3 sColor;
861 float afCubeValue[8];
862 float asEdgeVertex[12][3];
863 int pointIndex[12];
864
865 int fxScale = 1, fyScale = 1, fzScale = 1;
866 if ( cur_level != -1 )
867 {
868 fxScale = _emdata->get_xsize()/minvals[cur_level]->get_xsize();
869 fyScale = _emdata->get_ysize()/minvals[cur_level]->get_ysize();
870 fzScale = _emdata->get_zsize()/minvals[cur_level]->get_zsize();
871 for(iVertex = 0; iVertex < 8; iVertex++)
872 {
873 afCubeValue[iVertex] = _emdata->get_value_at( fxScale*(fX + a2fVertexOffset[iVertex][0]),
874 fyScale*(fY + a2fVertexOffset[iVertex][1]), fzScale*(fZ + a2fVertexOffset[iVertex][2]));
875 }
876 }
877 else
878 {
879 //Make a local copy of the values at the cube's corners
880 for(iVertex = 0; iVertex < 8; iVertex++)
881 {
882 afCubeValue[iVertex] = _emdata->get_value_at( fX + a2fVertexOffset[iVertex][0],
883 fY + a2fVertexOffset[iVertex][1], fZ + a2fVertexOffset[iVertex][2]);
884 }
885 }
886
887 //Find which vertices are inside of the surface and which are outside
888 iFlagIndex = 0;
889 for(iVertexTest = 0; iVertexTest < 8; iVertexTest++)
890 {
891 if (_surf_value >= 0 ){
892 if(afCubeValue[iVertexTest] <= _surf_value)
893 iFlagIndex |= 1<<iVertexTest;
894 }
895 else {
896 if(afCubeValue[iVertexTest] >= _surf_value)
897 iFlagIndex |= 1<<iVertexTest;
898 }
899 }
900
901 //Find which edges are intersected by the surface
902 iEdgeFlags = aiCubeEdgeFlags[iFlagIndex];
903
904 //If the cube is entirely inside or outside of the surface, then there will be no intersections
905 if(iEdgeFlags == 0) return;
906
907 //Find the point of intersection of the surface with each edge
908 //Then find the normal to the surface at those points
909 for(iEdge = 0; iEdge < 12; iEdge++)
910 {
911 //if there is an intersection on this edge
912 if(iEdgeFlags & (1<<iEdge))
913 {
914 fOffset = get_offset(afCubeValue[ a2iEdgeConnection[iEdge][0] ],
915 afCubeValue[ a2iEdgeConnection[iEdge][1] ], _surf_value);
916
917 if ( cur_level == -1 ){
918 asEdgeVertex[iEdge][0] = fX + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][0] + fOffset * a2fEdgeDirection[iEdge][0]) + 0.5f;
919 asEdgeVertex[iEdge][1] = fY + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][1] + fOffset * a2fEdgeDirection[iEdge][1]) + 0.5f;
920 asEdgeVertex[iEdge][2] = fZ + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][2] + fOffset * a2fEdgeDirection[iEdge][2]) + 0.5f;
921 } else {
922 asEdgeVertex[iEdge][0] = fxScale*(fX + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][0] + fOffset * a2fEdgeDirection[iEdge][0])) + 0.5f;
923 asEdgeVertex[iEdge][1] = fyScale*(fY + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][1] + fOffset * a2fEdgeDirection[iEdge][1])) + 0.5f;
924 asEdgeVertex[iEdge][2] = fzScale*(fZ + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][2] + fOffset * a2fEdgeDirection[iEdge][2])) + 0.5f;
925 }
926
927 pointIndex[iEdge] = get_edge_num(fX+edgeLookUp[iEdge][0], fY+edgeLookUp[iEdge][1], fZ+edgeLookUp[iEdge][2], edgeLookUp[iEdge][3]);
928 }
929 }
930
931 //Save voxel coords for later color processing
932 int vox[3] = {fX, fY, fZ};
933
934 //Draw the triangles that were found. There can be up to five per cube
935 for(iTriangle = 0; iTriangle < 5; iTriangle++)
936 {
937 if(a2iTriangleConnectionTable[iFlagIndex][3*iTriangle] < 0)
938 break;
939
940 float pts[3][3];
941 for(iCorner = 0; iCorner < 3; iCorner++)
942 {
943 iVertex = a2iTriangleConnectionTable[iFlagIndex][3*iTriangle+iCorner];
944 memcpy(&pts[iCorner][0], &asEdgeVertex[iVertex][0], 3*sizeof(float));
945 }
946
947
948
949 float v1[3] = {pts[1][0]-pts[0][0],pts[1][1]-pts[0][1],pts[1][2]-pts[0][2]};
950 float v2[3] = {pts[2][0]-pts[1][0],pts[2][1]-pts[1][1],pts[2][2]-pts[1][2]};
951
952 float n[3] = { v1[1]*v2[2] - v1[2]*v2[1], v1[2]*v2[0] - v1[0]*v2[2], v1[0]*v2[1] - v1[1]*v2[0] };
953
954
955 for(iCorner = 0; iCorner < 3; iCorner++)
956 {
957 // Without vertex normalization
958// iVertex = a2iTriangleConnectionTable[iFlagIndex][3*iTriangle+iCorner];
959// int ss = pp.elem();
960// pp.push_back_3(&pts[iCorner][0]);
961// nn.push_back_3(&n[0]);
962// ff.push_back(ss);
963
964// With vertex normalization
965 iVertex = a2iTriangleConnectionTable[iFlagIndex][3*iTriangle+iCorner];
966 map<int,int>::iterator it = point_map.find(pointIndex[iVertex]);
967 if ( it == point_map.end() ){
968 vv.push_back_3(&vox[0]);
969 int ss = pp.elem();
970 pp.push_back_3(&pts[iCorner][0]);
971 nn.push_back_3(&n[0]);
972 ff.push_back(ss);
973 point_map[pointIndex[iVertex]] = ss;
974 } else {
975 int idx = it->second;
976 ff.push_back(idx);
977 nn[idx] += n[0];
978 nn[idx+1] += n[1];
979 nn[idx+2] += n[2];
980 }
981 }
982 }
983}
984
985
986
987U3DWriter::U3DWriter() : DIFFUSE_COLOR_COUNT(1),
988 SPECULAR_COLOR_COUNT(1)
989{
990}
991U3DWriter::~U3DWriter() {}
992
993using std::ofstream;
994
995int U3DWriter::write(const string& filename) {
996 // Must open the ofstream in binary mode
997 ofstream of(filename.c_str(), ofstream::binary);
998 write(of);
999 of.close();
1000
1001 return 1;
1002}
1003
1004ostream& U3DWriter::write(ostream& os) {
1005 write_header(os);
1006 return os;
1007}
1008
1009unsigned int U3DWriter::size_of_in_bytes(){
1010 // this is just the size in bytes of all of the entries in this object that are written to the binary
1011 // output. This is based only on the behavior of write_header and needs to be udpated
1012 unsigned size = 4+2+2+4+4+8+4+8; // 36 bytes
1013 return size;
1014}
1015ostream& U3DWriter::write_header(ostream& os)
1016{
1017 // checks
1018 test_type_sizes();
1019
1020 U32 block_type_file_header = 0x00443355; // This is the block type tag of a file header, as taken form the ECMA spec
1021 write( os, block_type_file_header);
1022//
1023 I16 major_version = -1; // Compliance has not been tested for this encoder, so we must specify a value less than 0
1024 I16 minor_version = 0; // This is the version of this encoder, which we are calling '0'
1025 write( os, major_version);
1026 write( os, minor_version);
1027
1028 U32 profile_identifier = 0x00000000; // Base profile - no optional features ares used
1029 write( os, profile_identifier);
1030
1031 U32 declaration_size = size_of_in_bytes(); // This will have to be addressed at a later point, this is the size if the declaration block in bytes
1032 write( os, declaration_size);
1033
1034 U64 file_size = size_of_in_bytes(); // This is the size of the file in bytes
1035 write( os, file_size);
1036
1037 U32 character_encoding = 106; // UTF-8 MIB from the IANA
1038 write( os, character_encoding);
1039
1040 F64 unit_scaling = 1.0; // This should eventually be the correct scaling of the objects
1041 write( os, unit_scaling);
1042
1043
1044 return os;
1045}
1046
1047void U3DWriter::test_type_sizes()
1048{
1049 bool error = false;
1050 if (sizeof(F64) != 8 ){
1051 cout << "Error, size of double is not 64 bytes, it's " << sizeof(F64)*4 << endl;
1052 error = true;
1053 }
1054 if (sizeof(F32) != 4 ){
1055 cout << "Error, size of float is not 32 bytes, it's " << sizeof(F32)*4 << endl;
1056 error = true;
1057 }
1058 if (sizeof(U64) != 8) {
1059 cout << "Error, size of long unsigned int is not 64 bytes, it's " << sizeof(U64)*4 << endl;
1060 error = true;
1061 }
1062 if (sizeof(U32) != 4) {
1063 cout << "Error, size of unsigned int is not 32 bytes, it's " << sizeof(U32)*4 << endl;
1064 error = true;
1065 }
1066 if (sizeof(I16) != 2) {
1067 cout << "Error, size of short int is not 16 bytes, it's " << sizeof(I16)*4 << endl;
1068 error = true;
1069 }
1070 if (sizeof(U16) != 2) {
1071 cout << "Error, size of short unsigned int is not 16 bytes, it's " << sizeof(U16)*4 << endl;
1072 error = true;
1073 }
1074 if (sizeof(U8) != 1) {
1075 cout << "Error, size of unsigned char is not bytes, it's " << sizeof(U8)*4 << endl;
1076 error = true;
1077 }
1078
1079 if (error) {
1080 throw;
1081 }
1082}
1083
1084ostream& U3DWriter::write_clod_mesh_generator_node(ostream& os)
1085{
1086 /*
1087 CLOD Mesh Declaration
1088 */
1089 U32 block_type_clod_mesh_generator = 0xFFFFFF31; // This is the block type tag of the continuous level of detail mesh generator, as taken form the ECMA spec
1090 write( os, block_type_clod_mesh_generator);
1091
1092 string name("testmesh"); // The unique name, we get to make this up ourselves. It could be an empty string, we'd still have to call write_string(os,"")
1093 write(os,name);
1094
1095 U32 chain_index = 0x00000000; // the value of Chain Index shall be zero for this type - as in the SPEC
1096 write( os, chain_index);
1097
1098 /*
1099 Max Mesh Description
1100 */
1101 U32 mesh_attributes = 0x00000000; // Faces in the mesh have a normal index at each corner 0x00000001 is used to exlude normals POTENTIALLY USEFUL
1102 write(os,mesh_attributes);
1103 U32 face_count = ff.elem(); // The number of faces TO FILL IN LATER
1104 write(os,face_count);
1105 U32 position_count = pp.elem(); // The number of positions in the position array TO FILL IN LATER
1106 write(os,position_count);
1107 U32 normal_count = nn.elem(); // The number of normals in the normal array TO FILL IN LATER
1108 write(os,normal_count);
1109 U32 diffuse_color_count = DIFFUSE_COLOR_COUNT; // The number of colors in the diffuse color array TO FILL IN LATER
1110 write(os,diffuse_color_count);
1111 U32 specular_color_count = SPECULAR_COLOR_COUNT; // The number of colors in the specular color array TO FILL IN LATER
1112 write(os,specular_color_count);
1113 U32 texture_coord_count = 0x00000000; // The number of texture coordinates in teh texture coordinate array POTENTIALLY USEFUL
1114 write(os,texture_coord_count);
1115 U32 shading_count = 1; // The number of shading descriptions used in the mesh. This must correspond with the shader list in the shading group (see directly below we are using only one shader
1116 write(os,shading_count);
1117
1118 /*
1119 Shading Description
1120 */
1121 U32 shading_attributes = 0x00000003; // 0 means use neither diffuse or specular colors, 1 means use per vertex diffuse, 2 means use per vertex specular, 3 means use both POTENTIALLY USEFUL
1122 write(os,shading_attributes);
1123 U32 texture_layout_count = 0x00000000; // The number of texture layers used by this shader list
1124 write(os,texture_layout_count);
1125 U32 texture_coord_dimensions = 0x00000002; // The number of dimensions in the texture coordinate vector. I chose default to be 2. No particular reason. POTENTIALLY USEFUL
1126 write(os,texture_coord_dimensions);
1127 U32 original_shading_id = 0; // Original shading index for this shader list. Informative parameter only. This is shader 0
1128 write(os,original_shading_id);
1129
1130 /*
1131 CLOD Description - describes the range of resolutions available for the continuous level of detail mesh
1132 If there were more than one level of detail than these two numbers would be different
1133 */
1134 U32 minimum_resolution = position_count; // the number of positions in the base mesh
1135 write(os,minimum_resolution);
1136 U32 final_maximum_resolution = position_count; // the number of positions in the Max Mesh Description (by definition)
1137 write(os,final_maximum_resolution);
1138
1139 /*
1140 Resource Description
1141 */
1142
1143 /*
1144 Quality Factors
1145 for information only. Not used by the renderer. Helpful for conveying information to the user
1146 */
1147 U32 position_quality_factor = 0x00000000; // position quality factor. Descriptive information only
1148 write(os,position_quality_factor);
1149 U32 normal_quality_factor = 0x00000000; // normal quality factor. Descriptive information only
1150 write(os,normal_quality_factor);
1151 U32 texture_coord_quality_factor = 0x00000000; // texture coordinate quality factor. Descriptive information only
1152 write(os,texture_coord_quality_factor);
1153
1154 /*
1155 Inverse Quantization
1156 used to reconstruct floating point values that have been quantized.
1157 */
1158 F32 postion_inverse_quant = 1.0; // inverse quantization factor used in the reconstruction of the position vectors
1159 write(os,postion_inverse_quant);
1160 F32 normal_inverse_quant = 1.0; // inverse quantization factor used in the reconstruction of the normal vectors
1161 write(os,normal_inverse_quant);
1162 F32 texture_coord_inverse_quant = 1.0; // inverse quantization factor used in the reconstruction of the texture coordinates
1163 write(os,texture_coord_inverse_quant);
1164 F32 diffuse_color_inverse_quant = 1.0; // inverse quantization factor used in the reconstruction of the diffuse colors
1165 write(os,diffuse_color_inverse_quant);
1166 F32 specular_color_inverse_quant = 1.0; // inverse quantization factor used in the reconstruction of the specular colors
1167 write(os,specular_color_inverse_quant);
1168
1169 /*
1170 Resource Parameters
1171 parameters that help to define the conversion from the Author Mesh format to the Render Mesh format
1172 */
1173
1174 F32 normal_crease_parameter = 1.0; // A dot product value in the range -1 to 1, used to decide whether or not normals at the same position will be merged. 1.0 means never, -1.0 means always. 0 means the two normals must be separated by an angle no greater than 90 degrees to be merged. Think in angles.
1175 write(os,normal_crease_parameter);
1176 F32 normal_update_parameter = 0.0; // A strange a parameter that I couldn't make sense of - I think it means it will change the actual file itself if it encounters what it deems 'normal errors'
1177 write(os,normal_update_parameter);
1178 F32 normal_tolerance_parameter = 0.0; // Normals which are closer together than this value are considered equivalent in the Render Mesh. This is useful for compression
1179 write(os,normal_tolerance_parameter);
1180
1181 /*
1182 Skeleton description
1183 used in bones-based animation by the Animation Modifier
1184 */
1185 U32 bone_count = 0x00000000; // The number of bones associated with this mesh. We will always have 0, but this could change (unlikely).
1186 write(os,bone_count);
1187
1188 // WARNING - if bone_count is ever greater than one, then more writing would have to occur here
1189
1191 /*
1192 Base mesh is basically the minimum LOD mesh - it must exist
1193 */
1194 {
1195 U32 block_type_clod_base_mesh_continuation = 0xFFFFFF3B; // This is the block type tag of the CLOD Base Mesh Continuation, as taken form the ECMA spec
1196 write( os, block_type_clod_base_mesh_continuation);
1197
1198 write(os,name); // We use the same name as above
1199
1200 U32 chain_index = 0x00000000; // the value of Chain Index shall be zero for this type - as in the SPEC
1201 write( os, chain_index);
1202
1203 /*
1204 Base Mesh Description
1205 */
1206 U32 base_face_count = ff.elem(); // The number of faces in the base mesh TO FILL IN LATER
1207 write( os, base_face_count);
1208 U32 base_position_count = pp.elem(); // The number of positions used by base mesh in the position array TO FILL IN LATER
1209 write( os, base_position_count);
1210 U32 base_normal_count = nn.elem(); // The number of normals used by the base mesh in the normal array TO FILL IN LATER
1211 write( os, base_normal_count);
1212 U32 base_diffuse_color_count = DIFFUSE_COLOR_COUNT; // The number of diffuse colors used by the base mesh in the diffuse color array TO FILL IN LATER
1213 write( os, base_diffuse_color_count);
1214 U32 base_specular_color_count = SPECULAR_COLOR_COUNT; // The number of specular colors used by the base mesh in the specular color array TO FILL IN LATER
1215 write( os, base_specular_color_count);
1216 U32 base_texture_coord_count = 0x00000000; // The number of texture coordinates used by the base mesh in texture coordinate array TO FILL IN LATER
1217 write( os, base_texture_coord_count);
1218
1219 /*
1220 Base mesh data
1221 */
1222
1223 // Write position data
1224 F32* data = pp.get_data();
1225 for(unsigned int i = 0; i < pp.elem(); ++i) {
1226 write(os,data[i]);
1227 }
1228
1229 // Write normal data
1230 data = nn.get_data();
1231 for(unsigned int i = 0; i < nn.elem(); ++i) {
1232 write(os,data[i]);
1233 }
1234
1235 // Write Diffuse color, this is in rgba format
1236 F32 diffuse_rgba[4] = {1.0,0.0,0.0,1.0};
1237 for (unsigned int i = 0; i < 4; ++i) {
1238 write(os,diffuse_rgba[i]);
1239 }
1240
1241 // Write Specular color, this is in rgba format
1242 F32 specular_rgba[4] = {1.0,0.0,0.0,1.0};
1243 for (unsigned int i = 0; i < 4; ++i) {
1244 write(os,specular_rgba[i]);
1245 }
1246
1247 // THERE ARE NO TEXTURE COORDINATES, BUT IF THERE WERE WE WOULD HAVE TO WRITE THEM HERE
1248 // i.e. for i in range(0,base_texture_coord_count) write texture coords
1249
1250 // Write normal data
1251 U32* faces = ff.get_data();
1252 for(unsigned int i = 0; i < pp.elem(); i = i+3) {
1253 U32 shading_id = 0; // We only have one shader defined. This could could changed in future
1254 write(os,shading_id);
1255
1256 // write base corner info - there are always three corners
1257 for (unsigned int j =0; j < 3; ++j){
1258 U32 position_index = faces[i+j];
1259 write(os,position_index); // Write the position index
1260
1261 U32 normal_index = position_index;
1262 write(os,normal_index); // Write the normal index, it is exactly the same as the normal index!
1263
1264 U32 base_diffuse_color_index = 0; // only using one diffuse color
1265 write(os,base_diffuse_color_index);
1266
1267 U32 base_specular_color_index = 0; // only using one specular color
1268 write(os,base_specular_color_index);
1269
1270 // Would need to write texture data here if we were doing that
1271
1272 }
1273
1274 }
1275
1276 }
1277 return os;
1278}
1279
1280template<>
1281ostream& U3DWriter::write(ostream& os, const string& s )
1282{
1283 // WARNING - I AM NOT SURE IF THIS APPROACH WILL PRESENT UTF8 PROBLEMS.
1284 // See character_encoding in the file header writing module above
1285 test_type_sizes();
1286
1287 short unsigned int size = s.size(); // To write a string you must first place its
1288 write(os,size);
1289
1290 // Write the characters
1291 for(unsigned int i = 0; i<size; ++i) {
1292 write(os,s[i]);
1293 }
1294
1295 return os;
1296}
1297
1298#endif
EMAN::ColorRGBGenerator
Class to encapsulate an RGB color generator for marching cubes isosurface generator For now you can o...
Definition: marchingcubes.h:188
EMAN::ColorRGBGenerator::getRGBColor
float * getRGBColor(int x, int y, int z)
Generate a color based on pixel coords.
EMAN::ColorRGBGenerator::set_data
void set_data(EMData *data)
set the emdata
EMAN::ColorRGBGenerator::setNeedToRecolor
void setNeedToRecolor(bool recolor)
Definition: marchingcubes.h:248
EMAN::ColorRGBGenerator::set_cmap_data
void set_cmap_data(EMData *data)
Set min max data.
EMAN::ColorRGBGenerator::generateRadialColorMap
void generateRadialColorMap()
Generate a color map.
EMAN::ColorRGBGenerator::setMinMax
void setMinMax(float min, float max)
Set the mn max for cmap coloring.
Definition: marchingcubes.h:229
EMAN::Dict
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
EMAN::Dict::put
void put(const string &key, EMObject val)
Put the value/key pair into the dictionary probably better to just use operator[].
Definition: emobject.h:545
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
EMAN::Isosurface::_emdata
EMData * _emdata
Definition: isosurface.h:84
EMAN::MarchingCubes::marching_cube
void marching_cube(int fX, int fY, int fZ, const int cur_level)
Function for managing cases where a triangles can potentially be rendered Called by draw_cube.
EMAN::MarchingCubes::calculate_surface
void calculate_surface()
Calculate and generate the entire set of vertices and normals using current states Calls draw_cube(0,...
EMAN::MarchingCubes::set_surface_value
void set_surface_value(const float value)
Set Isosurface value.
EMAN::MarchingCubes::rgbgenerator
ColorRGBGenerator rgbgenerator
Color by radius generator.
Definition: marchingcubes.h:453
EMAN::MarchingCubes::cc
CustomVector< float > cc
Definition: marchingcubes.h:447
EMAN::MarchingCubes::get_normal
void get_normal(Vector3 &normal, int fX, int fY, int fZ)
Find the gradient of the scalar field at a point.
EMAN::MarchingCubes::calculate_min_max_vals
bool calculate_min_max_vals()
Calculate the min and max value trees Stores minimum and maximum cube neighborhood values in a tree s...
EMAN::MarchingCubes::get_offset
float get_offset(float fValue1, float fValue2, float fValueDesired)
Find the approximate point of intersection of the surface between two points with the values fValue1 ...
EMAN::MarchingCubes::set_data
void set_data(EMData *data)
Sets Voxel data for Isosurface implementation Calls calculate_min_max_vals which generates the tree o...
EMAN::MarchingCubes::maxvals
vector< EMData * > maxvals
Definition: marchingcubes.h:387
EMAN::MarchingCubes::draw_cube
void draw_cube(const int x, const int y, const int z, const int cur_level)
The main cube drawing function To start the process of generate triangles call with draw_cube(0,...
EMAN::MarchingCubes::drawing_level
int drawing_level
The "sampling rate".
Definition: marchingcubes.h:390
EMAN::MarchingCubes::nn
CustomVector< float > nn
Definition: marchingcubes.h:449
EMAN::MarchingCubes::get_isosurface
Dict get_isosurface()
Get the isosurface as dictionary Traverses the tree and marches the cubes.
EMAN::MarchingCubes::get_edge_num
int get_edge_num(int x, int y, int z, int edge)
Get edge num needs better commenting.
EMAN::MarchingCubes::ff
CustomVector< unsigned int > ff
Definition: marchingcubes.h:450
EMAN::MarchingCubes::MarchingCubes
MarchingCubes()
Default constructor.
EMAN::MarchingCubes::color_vertices
void color_vertices()
Color the vertices.
EMAN::MarchingCubes::~MarchingCubes
virtual ~MarchingCubes()
EMAN::MarchingCubes::clear_min_max_vals
void clear_min_max_vals()
Clear the minimum and maximum value search trees Frees memory in the minvals and maxvals.
EMAN::MarchingCubes::vv
CustomVector< int > vv
Definition: marchingcubes.h:448
EMAN::MarchingCubes::minvals
vector< EMData * > minvals
Vectors for storing the search trees.
Definition: marchingcubes.h:387
EMAN::MarchingCubes::pp
CustomVector< float > pp
.Custom vectors for storing points, normals and faces
Definition: marchingcubes.h:446
EMAN::MarchingCubes::point_map
map< int, int > point_map
Definition: marchingcubes.h:369
EMAN::Transform
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75
EMAN::U3DWriter::nn
CustomVector< F32 > nn
Definition: marchingcubes.h:611
EMAN::U3DWriter::ff
CustomVector< unsigned int > ff
Definition: marchingcubes.h:612
EMAN::U3DWriter::write_clod_mesh_generator_node
ostream & write_clod_mesh_generator_node(ostream &os)
EMAN::U3DWriter::test_type_sizes
void test_type_sizes()
EMAN::U3DWriter::write_header
ostream & write_header(ostream &)
EMAN::U3DWriter::U16
short unsigned int U16
Definition: marchingcubes.h:509
EMAN::U3DWriter::U64
long unsigned int U64
Definition: marchingcubes.h:505
EMAN::U3DWriter::pp
CustomVector< F32 > pp
Definition: marchingcubes.h:610
EMAN::U3DWriter::write
int write(const string &filename)
EMAN::U3DWriter::U8
unsigned char U8
Definition: marchingcubes.h:510
EMAN::U3DWriter::size_of_in_bytes
unsigned int size_of_in_bytes()
EMAN::U3DWriter::U32
unsigned int U32
Definition: marchingcubes.h:504
EMAN::Vector3::normalize
void normalize()
Definition: vecmath.h:276
sqrt
EMData * sqrt() const
return square root of current image
Definition: emdata_core.cpp:1084
get_ysize
int get_ysize() const
Get the image y-dimensional size.
Definition: emdata_metadata.h:553
get_zsize
int get_zsize() const
Get the image z-dimensional size.
Definition: emdata_metadata.h:562
get_xsize
int get_xsize() const
Get the image x-dimensional size.
Definition: emdata_metadata.h:544
set_data
void set_data(float *data, const int x, const int y, const int z)
Set the data explicitly data pointer must be allocated using malloc!
Definition: emdata_metadata.h:93
process
EMData * process(const string &processorname, const Dict &params=Dict()) const
Apply a processor with its parameters on a copy of this image, return result as a a new image.
NotExistingObjectException
#define NotExistingObjectException(objname, desc)
Definition: exception.h:130
ImageDimensionException
#define ImageDimensionException(desc)
Definition: exception.h:166
NullPointerException
#define NullPointerException(desc)
Definition: exception.h:241
EMAN
E2Exception class.
Definition: aligner.h:40
y
#define y(i, j)
Definition: projector.cpp:1516
x
#define x(i)
Definition: projector.cpp:1517
Generated by doxygen 1.9.2