NEURON
geometry3d.cpp
Go to the documentation of this file.
1 /*
2  This file contains code adapted from p3d.py in
3  http://code.google.com/p/pythonisosurfaces/source/checkout
4  which was released under the new BSD license.
5  accessed 31 July 2012
6 */
7 
8 #include <math.h>
9 #include <cstdio>
10 #include <cassert>
11 //#include <nrnpython.h>
12 
13  int geometry3d_find_triangles(double value0, double value1, double value2,
14  double value3, double value4, double value5,
15  double value6, double value7, double x0,
16  double x1, double y0, double y1, double z0,
17  double z1, double* out, int offset);
18 
19  double geometry3d_llgramarea(double* p0, double* p1, double* p2) ;
20  double geometry3d_sum_area_of_triangles(double* tri_vec, int len);
21 
22 inline double max(double a, double b) {
23  return a > b ? a : b;
24 }
25 
26 const int edgeTable[] = {
27 0x0, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
28 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
29 0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
30 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
31 0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
32 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
33 0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
34 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
35 0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
36 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
37 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
38 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
39 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
40 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
41 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
42 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
43 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
44 0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
45 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
46 0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
47 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
48 0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
49 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
50 0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
51 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
52 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
53 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
54 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
55 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
56 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
57 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
58 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0};
59 
60 /* CTNG:tritable */
61 const int triTable[256][16] = {
62 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
63 {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
64 {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
65 {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
66 {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
67 {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
68 {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
69 {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
70 {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
71 {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
72 {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
73 {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
74 {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
75 {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
76 {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
77 {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
78 {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
79 {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
80 {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
81 {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
82 {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
83 {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
84 {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
85 {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
86 {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
87 {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
88 {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
89 {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
90 {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
91 {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
92 {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
93 {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
94 {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
95 {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
96 {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
97 {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
98 {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
99 {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
100 {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
101 {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
102 {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
103 {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
104 {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
105 {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
106 {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
107 {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
108 {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
109 {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
110 {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
111 {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
112 {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
113 {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
114 {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
115 {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
116 {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
117 {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
118 {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
119 {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
120 {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
121 {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
122 {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
123 {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
124 {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
125 {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
126 {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
127 {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
128 {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
129 {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
130 {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
131 {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
132 {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
133 {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
134 {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
135 {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
136 {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
137 {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
138 {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
139 {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
140 {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
141 {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
142 {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
143 {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
144 {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
145 {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
146 {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
147 {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
148 {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
149 {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
150 {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
151 {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
152 {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
153 {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
154 {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
155 {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
156 {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
157 {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
158 {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
159 {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
160 {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
161 {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
162 {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
163 {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
164 {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
165 {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
166 {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
167 {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
168 {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
169 {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
170 {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
171 {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
172 {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
173 {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
174 {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
175 {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
176 {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
177 {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
178 {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
179 {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
180 {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
181 {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
182 {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
183 {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
184 {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
185 {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
186 {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
187 {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
188 {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
189 {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
190 {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
191 {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
192 {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
193 {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
194 {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
195 {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
196 {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
197 {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
198 {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
199 {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
200 {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
201 {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
202 {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
203 {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
204 {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
205 {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
206 {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
207 {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
208 {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
209 {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
210 {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
211 {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
212 {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
213 {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
214 {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
215 {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
216 {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
217 {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
218 {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
219 {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
220 {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
221 {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
222 {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
223 {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
224 {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
225 {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
226 {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
227 {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
228 {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
229 {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
230 {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
231 {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
232 {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
233 {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
234 {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
235 {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
236 {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
237 {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
238 {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
239 {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
240 {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
241 {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
242 {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
243 {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
244 {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
245 {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
246 {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
247 {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
248 {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
249 {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
250 {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
251 {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
252 {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
253 {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
254 {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
255 {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
256 {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
257 {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
258 {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
259 {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
260 {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
261 {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
262 {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
263 {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
264 {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
265 {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
266 {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
267 {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
268 {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
269 {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
270 {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
271 {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
272 {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
273 {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
274 {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
275 {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
276 {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
277 {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
278 {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
279 {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
280 {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
281 {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
282 {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
283 {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
284 {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
285 {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
286 {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
287 {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
288 {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
289 {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
290 {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
291 {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
292 {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
293 {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
294 {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
295 {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
296 {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
297 {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
298 {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
299 {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
300 {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
301 {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
302 {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
303 {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
304 {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
305 {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
306 {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
307 {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
308 {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
309 {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
310 {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
311 {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
312 {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
313 {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
314 {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
315 {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
316 {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
317 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}};
318 
319 
320 
321 /* CTNG:interpedge */
322 void geometry3d_vi(double* p1, double* p2, double v1, double v2, double* out) {
323  if (fabs(v1) < 0.000000000001) {
324  out[0] = p1[0]; out[1] = p1[1]; out[2] = p1[2];
325  return;
326  }
327  if (fabs(v2) < 0.000000000001) {
328  out[0] = p2[0]; out[1] = p2[1]; out[2] = p2[2];
329  return;
330  }
331  double delta_v = v1 - v2;
332  if (fabs(delta_v) < 0.0000000001) {
333  out[0] = p1[0]; out[1] = p1[1]; out[2] = p1[2];
334  return;
335  }
336  double mu = v1 / delta_v;
337  if (isnan(mu)) {
338  printf("geometry3d_vi error. delta_v = %g, v1 = %g, v2 = %g\n", delta_v, v1, v2);
339  }
340  out[0] = p1[0] + mu * (p2[0] - p1[0]);
341  out[1] = p1[1] + mu * (p2[1] - p1[1]);
342  out[2] = p1[2] + mu * (p2[2] - p1[2]);
343 }
344 
345 int geometry3d_find_triangles(double value0, double value1, double value2, double value3, double value4, double value5, double value6, double value7, double x0, double x1, double y0, double y1, double z0, double z1, double* out, int offset) {
346  out = out + offset;
347  double position[8][3] = {{x0, y0, z0},
348  {x1, y0, z0},
349  {x1, y1, z0},
350  {x0, y1, z0},
351  {x0, y0, z1},
352  {x1, y0, z1},
353  {x1, y1, z1},
354  {x0, y1, z1}};
355 
356  /* CTNG:domarch */
357  int cubeIndex = 0;
358  if (value0 < 0) cubeIndex |= 1;
359  if (value1 < 0) cubeIndex |= 2;
360  if (value2 < 0) cubeIndex |= 4;
361  if (value3 < 0) cubeIndex |= 8;
362  if (value4 < 0) cubeIndex |= 16;
363  if (value5 < 0) cubeIndex |= 32;
364  if (value6 < 0) cubeIndex |= 64;
365  if (value7 < 0) cubeIndex |= 128;
366 
367  int et = edgeTable[cubeIndex];
368 
369  if (et == 0) return 0;
370 
371  double vertexList[12][3];
372  if (et & 1) geometry3d_vi(position[0], position[1], value0, value1, vertexList[0]);
373  if (et & 2) geometry3d_vi(position[1], position[2], value1, value2, vertexList[1]);
374  if (et & 4) geometry3d_vi(position[2], position[3], value2, value3, vertexList[2]);
375  if (et & 8) geometry3d_vi(position[3], position[0], value3, value0, vertexList[3]);
376  if (et & 16) geometry3d_vi(position[4], position[5], value4, value5, vertexList[4]);
377  if (et & 32) geometry3d_vi(position[5], position[6], value5, value6, vertexList[5]);
378  if (et & 64) geometry3d_vi(position[6], position[7], value6, value7, vertexList[6]);
379  if (et & 128) geometry3d_vi(position[7], position[4], value7, value4, vertexList[7]);
380  if (et & 256) geometry3d_vi(position[0], position[4], value0, value4, vertexList[8]);
381  if (et & 512) geometry3d_vi(position[1], position[5], value1, value5, vertexList[9]);
382  if (et & 1024) geometry3d_vi(position[2], position[6], value2, value6, vertexList[10]);
383  if (et & 2048) geometry3d_vi(position[3], position[7], value3, value7, vertexList[11]);
384 
385  int const * const tt = triTable[cubeIndex];
386  int i, j, k, count;
387  for (i = 0, count = 0; i < 16; i += 3, count++) {
388  if (tt[i] == -1) break;
389  for (k = 0; k < 3; k++) {
390  for (j = 0; j < 3; j++) out[j] = vertexList[tt[i + k]][j];
391  out += 3;
392  }
393  }
394  return count;
395 }
396 
397 
398 
399 double geometry3d_llgramarea(double* p0, double* p1, double* p2) {
400  /* setup the vectors */
401  double a[] = {p0[0] - p1[0], p0[1] - p1[1], p0[2] - p1[2]};
402  double b[] = {p0[0] - p2[0], p0[1] - p2[1], p0[2] - p2[2]};
403 
404  /* take the cross-product */
405  double cpx = a[1] * b[2] - a[2] * b[1];
406  double cpy = a[2] * b[0] - a[0] * b[2];
407  double cpz = a[0] * b[1] - a[1] * b[0];
408  return sqrt(cpx * cpx + cpy * cpy + cpz * cpz);
409 }
410 
411 
412 
413 double geometry3d_sum_area_of_triangles(double* tri_vec, int len) {
414  double sum = 0;
415  for (int i = 0; i < len; i += 9) {
416  sum += geometry3d_llgramarea(tri_vec + i, tri_vec + i + 3, tri_vec + i + 6);
417  }
418  return sum / 2.;
419 }
420 
421 
422 
423 /*****************************************************************************
424  * this contains cone, and cylinder code translated and modified
425  * from Barbier and Galin 2004's example code
426  * see /u/ramcd/spatial/experiments/one_time_tests/2012-06-28/cone.cpp
427  * on 7 june 2012, the original code was available online at
428  * http://jgt.akpeters.com/papers/BarbierGalin04/Cone-Sphere.zip
429  ****************************************************************************/
430 
431 
432 class geometry3d_Cylinder {
433  public:
434  geometry3d_Cylinder(double x0, double y0, double z0, double x1, double y1, double z1, double r);
435  //~geometry3d_Cylinder();
436  double signed_distance(double px, double py, double pz);
437  private:
438  //double x0, y0, z0, x1, y1, z1, r;
439  double r, rr, axisx, axisy, axisz, cx, cy, cz, h;
440 };
441 
442 geometry3d_Cylinder::geometry3d_Cylinder(double x0, double y0, double z0,
443  double x1, double y1, double z1, double r) :
444  r(r), cx((x0 + x1) / 2.), cy((y0 + y1) / 2.), cz((z0 + z1) / 2.),
445  rr(r * r) {
446  axisx = x1 - x0;
447  axisy = y1 - y0;
448  axisz = z1 - z0;
449  double axislength = sqrt(axisx * axisx + axisy * axisy + axisz * axisz);
450  axisx /= axislength;
451  axisy /= axislength;
452  axisz /= axislength;
453  h = axislength / 2.;
454 }
455 
456 double geometry3d_Cylinder::signed_distance(double px, double py, double pz) {
457  double nx, ny, nz, y, yy, xx, x;
458  nx = px - cx;
459  ny = py - cy;
460  nz = pz - cz;
461  y = fabs(axisx * nx + axisy * ny + axisz * nz);
462  yy = y * y;
463  xx = nx * nx + ny * ny + nz * nz - yy;
464  if (y < h) {
465  return max(-fabs(h - y), sqrt(xx) - r);
466  } else {
467  y -= h;
468  yy = y * y;
469  if (xx < rr) {
470  return fabs(y);
471  } else {
472  x = sqrt(xx) - r;
473  return sqrt(yy + x * x);
474  }
475  }
476 }
477 
478  void* geometry3d_new_Cylinder(double x0, double y0, double z0, double x1, double y1, double z1, double r) {
479  return new geometry3d_Cylinder(x0, y0, z0, x1, y1, z1, r);
480  }
481  void geometry3d_delete_Cylinder(void* ptr) {
482  delete (geometry3d_Cylinder*) ptr;
483  }
484  // TODO: add validation for ptr
485  double geometry3d_Cylinder_signed_distance(void* ptr, double px, double py, double pz) {
486  return ((geometry3d_Cylinder*) ptr) -> signed_distance(px, py, pz);
487  }
488 
489 
490 
491 
492 class geometry3d_Cone {
493  public:
494  geometry3d_Cone(double x0, double y0, double z0, double r0, double x1, double y1, double z1, double r1);
495  double signed_distance(double px, double py, double pz);
496  private:
497  double axisx, axisy, axisz, cx, cy, cz, h, rra, rrb, conelength;
498  double side1, side2, x0, y0, z0, r0, axislength;
499 };
500 
501 geometry3d_Cone::geometry3d_Cone(double x0, double y0, double z0, double r0,
502  double x1, double y1, double z1, double r1) :
503  cx((x0 + x1) / 2.), cy((y0 + y1) / 2.), cz((z0 + z1) / 2.),
504  rra(r0 * r0), rrb(r1 * r1), x0(x0), y0(y0), z0(z0), r0(r0) {
505  // TODO: these are preconditions; the python assures them, but could/should
506  // take care of that here
507  assert(r1 <= r0);
508  assert(r1 >= 0);
509  axisx = x1 - x0;
510  axisy = y1 - y0;
511  axisz = z1 - z0;
512  axislength = sqrt(axisx * axisx + axisy * axisy + axisz * axisz);
513  axisx /= axislength;
514  axisy /= axislength;
515  axisz /= axislength;
516  h = axislength / 2.;
517  rra = r0 * r0;
518  rrb = r1 * r1;
519  conelength = sqrt((r1 - r0) * (r1 - r0) + axislength * axislength);
520  side1 = (r1 - r0) / conelength;
521  side2 = axislength / conelength;
522 }
523 
524 double geometry3d_Cone::signed_distance(double px, double py, double pz) {
525  double nx, ny, nz, y, yy, xx, x, rx, ry;
526  nx = px - x0;
527  ny = py - y0;
528  nz = pz - z0;
529  y = axisx * nx + axisy * ny + axisz * nz;
530  yy = y * y;
531  xx = nx * nx + ny * ny + nz * nz - yy;
532  // in principle, xx >= 0, but roundoff errors may cause trouble
533  if (xx < 0) xx = 0;
534  if (y < 0) {
535  if (xx < rra) return -y;
536  x = sqrt(xx) - r0;
537  return sqrt(x * x + yy);
538  } else if (xx < rrb) {
539  return y - axislength;
540  } else {
541  x = sqrt(xx) - r0;
542  if (y < 0) {
543  if (x < 0) return y;
544  return sqrt(x * x + yy);
545  } else {
546  ry = x * side1 + y * side2;
547  if (ry < 0) return sqrt(x * x + yy);
548  rx = x * side2 - y * side1;
549  if (ry > conelength) {
550  ry -= conelength;
551  return sqrt(rx * rx + ry * ry);
552  } else {
553  return rx;
554  }
555  }
556  }
557 }
558 
559 
560  void* geometry3d_new_Cone(double x0, double y0, double z0, double r0, double x1, double y1, double z1, double r1) {
561  return new geometry3d_Cone(x0, y0, z0, r0, x1, y1, z1, r1);
562  }
563  void geometry3d_delete_Cone(void* ptr) {
564  delete (geometry3d_Cone*) ptr;
565  }
566  // TODO: add validation for ptr
567  double geometry3d_Cone_signed_distance(void* ptr, double px, double py, double pz) {
568  return ((geometry3d_Cone*) ptr) -> signed_distance(px, py, pz);
569  }
570 
571 
572 class geometry3d_Sphere {
573  public:
574  geometry3d_Sphere(double x, double y, double z, double r);
575  double signed_distance(double px, double py, double pz);
576  private:
577  double x, y, z, r;
578 };
579 
580 geometry3d_Sphere::geometry3d_Sphere(double x, double y, double z, double r) :
581  x(x), y(y), z(z), r(r) {}
582 
583 double geometry3d_Sphere::signed_distance(double px, double py, double pz) {
584  return sqrt(pow(x - px, 2) + pow(y - py, 2) + pow(z - pz, 2)) - r;
585 }
586 
587  void* geometry3d_new_Sphere(double x, double y, double z, double r) {
588  return new geometry3d_Sphere(x, y, z, r);
589  }
590  void geometry3d_delete_Sphere(void* ptr) {
591  delete (geometry3d_Sphere*) ptr;
592  }
593  // TODO: add validation for ptr
594  double geometry3d_Sphere_signed_distance(void* ptr, double px, double py, double pz) {
595  return ((geometry3d_Sphere*) ptr) -> signed_distance(px, py, pz);
596  }
597 
598 class geometry3d_Plane {
599  public:
600  geometry3d_Plane(double x, double y, double z, double nx, double ny, double nz);
601  double signed_distance(double px, double py, double pz);
602  private:
603  double nx, ny, nz;
604  double d, mul;
605 };
606 
607 geometry3d_Plane::geometry3d_Plane(double x, double y, double z, double nx, double ny, double nz) :
608  nx(nx), ny(ny), nz(nz), d(-(nx * x + ny * y + nz * z)),
609  mul(1. / sqrt(nx * nx + ny * ny + nz * nz)) {}
610 
611 double geometry3d_Plane::signed_distance(double px, double py, double pz) {
612  return (nx * px + ny * py + nz * pz + d) * mul;
613 }
614 
615  void* geometry3d_new_Plane(double x, double y, double z, double nx, double ny, double nz) {
616  return new geometry3d_Plane(x, y, z, nx, ny, nz);
617  }
618  void geometry3d_delete_Plane(void* ptr) {
619  delete (geometry3d_Plane*) ptr;
620  }
621  // TODO: add validation for ptr
622  double geometry3d_Plane_signed_distance(void* ptr, double px, double py, double pz) {
623  return ((geometry3d_Plane*) ptr) -> signed_distance(px, py, pz);
624  }
625 
626 /*
627  PyObject* nrnpy_pyCallObject(PyObject*, PyObject*);
628 
629  void print_numbers(PyObject *p) {
630  for (Py_ssize_t i = 0; i< PyList_Size(p); i++) {
631  PyObject* obj = PyList_GetItem(p, i);
632  printf("%g ", PyFloat_AsDouble(obj));
633  }
634  printf("\n");
635  }
636 
637  // TODO: it would be nice to remove the python dependence, because that
638  // limits us to mostly single threaded due to the global interpreter
639  // lock
640  int geometry3d_process_cell(int i, int j, int k, PyObject* objects, double* xs, double* ys, double* zs, double* tridata, int start) {
641  double x, y, z, x1, y1, z1, xx, yy, zz;
642  x = xs[i]; y = ys[j]; z = zs[k];
643  x1 = xs[i + 1]; y1 = ys[j + 1]; z1 = zs[k + 1];
644  double value[8], current_value;
645  PyObject* result;
646  PyObject* obj;
647  printf("inside process_cell\n");
648 
649  // march around the cube
650  for (int m = 0; m < 8; m++) {
651  // NOTE: only describing changes from previous case
652  switch(m) {
653  case 0: xx = x; yy = y; zz = z; break;
654  case 1: xx = x1; break;
655  case 2: yy = y1; break;
656  case 3: xx = x; break;
657  case 4: yy = y; zz = z1; break;
658  case 5: xx = x1; break;
659  case 6: yy = y1; break;
660  case 7: xx = x; break;
661  }
662  printf("phase 0, len(objects) = %ld\n", PyList_Size(objects));
663  obj = PyList_GetItem(objects, 0);
664  printf("phase 0a, obj = %x\n", (void*) obj);
665  result = PyEval_CallMethod(obj, "distance", "ddd", xx, yy, zz);
666  //Py_DECREF(obj);
667  printf("phase 1\n");
668  current_value = PyFloat_AsDouble(result);
669  //Py_DECREF(result);
670  for (Py_ssize_t n = 1; n < PyList_Size(objects); n++) {
671  printf("phase 2, n = %ld\n", n);
672  obj = PyList_GetItem(objects, n);
673  result = PyEval_CallMethod(obj, "distance", "ddd", xx, yy, zz);
674  Py_DECREF(obj);
675  current_value = min(current_value, PyFloat_AsDouble(result));
676  //Py_DECREF(result);
677  }
678  value[m] = current_value;
679  }
680  printf("finishing up; start = %d\n", start);
681  return start + 9 * geometry3d_find_triangles(value[0], value[1], value[2], value[3], value[4], value[5], value[6], value[7], x, x1, y, y1, z, z1, tridata, start);
682  }
683 
684 
685  int geometry3d_test_call_function(PyObject* obj) {
686  printf("inside\n");
687  if (obj == NULL) printf("obj is NULL\n");
688  Py_INCREF(obj);
689  PyEval_CallObject(obj, obj);
690  return 0;
691  }
692 
693  int geometry3d_test_call_function3(int (*func) (void)) {
694  return func();
695  }
696 
697  int geometry3d_test_call_method(PyObject* list, PyObject* obj) {
698  PyEval_CallMethod(list, "insert", "O", obj);
699  return 0;
700  }
701 
702  // this works, but is slower than the python version
703  int geometry3d_contains_surface(int i, int j, int k, double (*objdist)(double, double, double), double* xs, double* ys, double* zs, double dx, double r_inner, double r_outer) {
704  bool has_neg = false, has_pos = false;
705  double xbar = xs[i] + dx / 2;
706  double ybar = ys[j] + dx / 2;
707  double zbar = zs[k] + dx / 2;
708 
709  double d = fabs(objdist(xbar, ybar, zbar));
710  //if (i == 586 && j == 2169 && k == 83) {printf("at magic point, d=%g\n", d);}
711  //printf("sphere test: d = %g, r_inner = %g, r_outer = %g\n", d, r_inner, r_outer);
712  if (d <= r_inner) return 1;
713  if (d >= r_outer) return 0;
714 // // spheres alone are indeterminant; check corners
715 // for (int di = 0; di < 2; di++) {
716 // for (int dj = 0; dj < 2; dj++) {
717 // for (int dk = 0; dk < 2; dk++) {
718 // d = objdist(xs[i + di], ys[j + dj], zs[k + dk]);
719 // //printf("d = %g\n", d);
720 // if (d <= 0) has_neg = true;
721 // if (d >= 0) has_pos = true;
722 // if (has_neg && has_pos) return 1;
723 // }
724 // }
725 // }
726 
727  // spheres alone are indeterminant; check corners
728  for (double* x = xs + i; x < xs + i + 2; x++) {
729  for (double* y = ys + j; y < ys + j + 2; y++) {
730  for (double* z = zs + k; z < zs + k + 2; z++) {
731  d = objdist(*x, *y, *z);
732  if (d <= 0) has_neg = true;
733  if (d >= 0) has_pos = true;
734  if (has_neg && has_pos) return 1;
735  }
736  }
737  }
738 
739  return 0;
740  }
741 
742 }
743 */
geometry3d_Cylinder::rr
double rr
Definition: geometry3d.cpp:439
geometry3d_Sphere
Definition: geometry3d.cpp:572
max
double max(double a, double b)
Definition: geometry3d.cpp:22
assert
#define assert(ex)
Definition: hocassrt.h:26
geometry3d_Cone::geometry3d_Cone
geometry3d_Cone(double x0, double y0, double z0, double r0, double x1, double y1, double z1, double r1)
Definition: geometry3d.cpp:501
geometry3d_new_Plane
void * geometry3d_new_Plane(double x, double y, double z, double nx, double ny, double nz)
Definition: geometry3d.cpp:615
geometry3d_vi
void geometry3d_vi(double *p1, double *p2, double v1, double v2, double *out)
Definition: geometry3d.cpp:322
geometry3d_Cylinder::cx
double cx
Definition: geometry3d.cpp:439
geometry3d_Cylinder::geometry3d_Cylinder
geometry3d_Cylinder(double x0, double y0, double z0, double x1, double y1, double z1, double r)
Definition: geometry3d.cpp:442
geometry3d_Plane
Definition: geometry3d.cpp:598
geometry3d_find_triangles
int geometry3d_find_triangles(double value0, double value1, double value2, double value3, double value4, double value5, double value6, double value7, double x0, double x1, double y0, double y1, double z0, double z1, double *out, int offset)
Definition: geometry3d.cpp:345
geometry3d_Cone::side2
double side2
Definition: geometry3d.cpp:498
geometry3d_Cone::side1
double side1
Definition: geometry3d.cpp:498
geometry3d_Plane::nx
double nx
Definition: geometry3d.cpp:603
geometry3d_Sphere_signed_distance
double geometry3d_Sphere_signed_distance(void *ptr, double px, double py, double pz)
Definition: geometry3d.cpp:594
geometry3d_Cone
Definition: geometry3d.cpp:492
geometry3d_new_Cone
void * geometry3d_new_Cone(double x0, double y0, double z0, double r0, double x1, double y1, double z1, double r1)
Definition: geometry3d.cpp:560
geometry3d_llgramarea
double geometry3d_llgramarea(double *p0, double *p1, double *p2)
Definition: geometry3d.cpp:399
geometry3d_Cone::z0
double z0
Definition: geometry3d.cpp:498
geometry3d_Sphere::y
double y
Definition: geometry3d.cpp:577
geometry3d_Cylinder_signed_distance
double geometry3d_Cylinder_signed_distance(void *ptr, double px, double py, double pz)
Definition: geometry3d.cpp:485
geometry3d_delete_Cone
void geometry3d_delete_Cone(void *ptr)
Definition: geometry3d.cpp:563
k
static philox4x32_key_t k
Definition: nrnran123.cpp:11
geometry3d_Cone::h
double h
Definition: geometry3d.cpp:497
geometry3d_delete_Cylinder
void geometry3d_delete_Cylinder(void *ptr)
Definition: geometry3d.cpp:481
edgeTable
const int edgeTable[]
Definition: geometry3d.cpp:26
geometry3d_Plane_signed_distance
double geometry3d_Plane_signed_distance(void *ptr, double px, double py, double pz)
Definition: geometry3d.cpp:622
mul
#define mul
Definition: redef.h:93
geometry3d_Sphere::x
double x
Definition: geometry3d.cpp:577
geometry3d_new_Sphere
void * geometry3d_new_Sphere(double x, double y, double z, double r)
Definition: geometry3d.cpp:587
geometry3d_Cylinder::r
double r
Definition: geometry3d.cpp:439
geometry3d_Cylinder::h
double h
Definition: geometry3d.cpp:439
geometry3d_Plane::mul
double mul
Definition: geometry3d.cpp:604
printf
#define printf
Definition: mwprefix.h:26
geometry3d_Cone::rrb
double rrb
Definition: geometry3d.cpp:497
triTable
const int triTable[256][16]
Definition: geometry3d.cpp:61
geometry3d_Plane::ny
double ny
Definition: geometry3d.cpp:603
geometry3d_Cylinder::cz
double cz
Definition: geometry3d.cpp:439
geometry3d_Sphere::r
double r
Definition: geometry3d.cpp:577
geometry3d_delete_Sphere
void geometry3d_delete_Sphere(void *ptr)
Definition: geometry3d.cpp:590
pow
pow
Definition: extdef.h:3
j
size_t j
Definition: nrngsl_real_radix2.cpp:56
geometry3d_Cone::axislength
double axislength
Definition: geometry3d.cpp:498
geometry3d_Cone::axisy
double axisy
Definition: geometry3d.cpp:497
geometry3d_Plane::d
double d
Definition: geometry3d.cpp:604
geometry3d_Cone::axisz
double axisz
Definition: geometry3d.cpp:497
geometry3d_delete_Plane
void geometry3d_delete_Plane(void *ptr)
Definition: geometry3d.cpp:618
geometry3d_Cylinder::cy
double cy
Definition: geometry3d.cpp:439
geometry3d_Cone::y0
double y0
Definition: geometry3d.cpp:498
geometry3d_Cone::r0
double r0
Definition: geometry3d.cpp:498
geometry3d_Plane::nz
double nz
Definition: geometry3d.cpp:603
geometry3d_Cylinder::axisz
double axisz
Definition: geometry3d.cpp:439
geometry3d_Sphere::geometry3d_Sphere
geometry3d_Sphere(double x, double y, double z, double r)
Definition: geometry3d.cpp:580
sqrt
sqrt
Definition: extdef.h:3
geometry3d_Cylinder
Definition: geometry3d.cpp:432
geometry3d_Cone::x0
double x0
Definition: geometry3d.cpp:498
geometry3d_Cylinder::axisx
double axisx
Definition: geometry3d.cpp:439
geometry3d_sum_area_of_triangles
double geometry3d_sum_area_of_triangles(double *tri_vec, int len)
Definition: geometry3d.cpp:413
geometry3d_Cone::signed_distance
double signed_distance(double px, double py, double pz)
Definition: geometry3d.cpp:524
i
#define i
Definition: md1redef.h:12
geometry3d_Cone::rra
double rra
Definition: geometry3d.cpp:497
geometry3d_Cylinder::axisy
double axisy
Definition: geometry3d.cpp:439
geometry3d_Cone_signed_distance
double geometry3d_Cone_signed_distance(void *ptr, double px, double py, double pz)
Definition: geometry3d.cpp:567
geometry3d_new_Cylinder
void * geometry3d_new_Cylinder(double x0, double y0, double z0, double x1, double y1, double z1, double r)
Definition: geometry3d.cpp:478
fabs
fabs
Definition: extdef.h:3
geometry3d_Plane::geometry3d_Plane
geometry3d_Plane(double x, double y, double z, double nx, double ny, double nz)
Definition: geometry3d.cpp:607
geometry3d_Plane::signed_distance
double signed_distance(double px, double py, double pz)
Definition: geometry3d.cpp:611
geometry3d_Sphere::signed_distance
double signed_distance(double px, double py, double pz)
Definition: geometry3d.cpp:583
geometry3d_Cone::axisx
double axisx
Definition: geometry3d.cpp:497
geometry3d_Sphere::z
double z
Definition: geometry3d.cpp:577
geometry3d_Cone::conelength
double conelength
Definition: geometry3d.cpp:497
geometry3d_Cylinder::signed_distance
double signed_distance(double px, double py, double pz)
Definition: geometry3d.cpp:456