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