ClangFormat: apply to source, most of intern
[blender.git] / intern / dualcon / intern / Projections.cpp
1 /*
2  * This program is free software; you can redistribute it and/or
3  * modify it under the terms of the GNU General Public License
4  * as published by the Free Software Foundation; either version 2
5  * of the License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program; if not, write to the Free Software Foundation,
14  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15  */
16
17 #ifdef WITH_CXX_GUARDEDALLOC
18 #  include "MEM_guardedalloc.h"
19 #endif
20
21 #include <math.h>
22 #include "Projections.h"
23
24 const int vertmap[8][3] = {
25     {0, 0, 0},
26     {0, 0, 1},
27     {0, 1, 0},
28     {0, 1, 1},
29     {1, 0, 0},
30     {1, 0, 1},
31     {1, 1, 0},
32     {1, 1, 1},
33 };
34
35 const int centmap[3][3][3][2] = {
36     {{{0, 0}, {0, 1}, {1, 1}}, {{0, 2}, {0, 3}, {1, 3}}, {{2, 2}, {2, 3}, {3, 3}}},
37
38     {{{0, 4}, {0, 5}, {1, 5}}, {{0, 6}, {0, 7}, {1, 7}}, {{2, 6}, {2, 7}, {3, 7}}},
39
40     {{{4, 4}, {4, 5}, {5, 5}}, {{4, 6}, {4, 7}, {5, 7}}, {{6, 6}, {6, 7}, {7, 7}}}};
41
42 const int edgemap[12][2] = {
43     {0, 4},
44     {1, 5},
45     {2, 6},
46     {3, 7},
47     {0, 2},
48     {1, 3},
49     {4, 6},
50     {5, 7},
51     {0, 1},
52     {2, 3},
53     {4, 5},
54     {6, 7},
55 };
56
57 const int facemap[6][4] = {
58     {0, 1, 2, 3},
59     {4, 5, 6, 7},
60     {0, 1, 4, 5},
61     {2, 3, 6, 7},
62     {0, 2, 4, 6},
63     {1, 3, 5, 7},
64 };
65
66 /**
67  * Method to perform cross-product
68  */
69 static void crossProduct(int64_t res[3], const int64_t a[3], const int64_t b[3])
70 {
71   res[0] = a[1] * b[2] - a[2] * b[1];
72   res[1] = a[2] * b[0] - a[0] * b[2];
73   res[2] = a[0] * b[1] - a[1] * b[0];
74 }
75
76 static void crossProduct(double res[3], const double a[3], const double b[3])
77 {
78   res[0] = a[1] * b[2] - a[2] * b[1];
79   res[1] = a[2] * b[0] - a[0] * b[2];
80   res[2] = a[0] * b[1] - a[1] * b[0];
81 }
82
83 /**
84  * Method to perform dot product
85  */
86 static int64_t dotProduct(const int64_t a[3], const int64_t b[3])
87 {
88   return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
89 }
90
91 static void normalize(double a[3])
92 {
93   double mag = a[0] * a[0] + a[1] * a[1] + a[2] * a[2];
94   if (mag > 0) {
95     mag = sqrt(mag);
96     a[0] /= mag;
97     a[1] /= mag;
98     a[2] /= mag;
99   }
100 }
101
102 /* Create projection axes for cube+triangle intersection testing.
103  *    0, 1, 2: cube face normals
104  *
105  *          3: triangle normal
106  *
107  *    4, 5, 6,
108  *    7, 8, 9,
109  * 10, 11, 12: cross of each triangle edge vector with each cube
110  *             face normal
111  */
112 static void create_projection_axes(int64_t axes[NUM_AXES][3], const int64_t tri[3][3])
113 {
114   /* Cube face normals */
115   axes[0][0] = 1;
116   axes[0][1] = 0;
117   axes[0][2] = 0;
118   axes[1][0] = 0;
119   axes[1][1] = 1;
120   axes[1][2] = 0;
121   axes[2][0] = 0;
122   axes[2][1] = 0;
123   axes[2][2] = 1;
124
125   /* Get triangle edge vectors */
126   int64_t tri_edges[3][3];
127   for (int i = 0; i < 3; i++) {
128     for (int j = 0; j < 3; j++)
129       tri_edges[i][j] = tri[(i + 1) % 3][j] - tri[i][j];
130   }
131
132   /* Triangle normal */
133   crossProduct(axes[3], tri_edges[0], tri_edges[1]);
134
135   // Face edges and triangle edges
136   int ct = 4;
137   for (int i = 0; i < 3; i++) {
138     for (int j = 0; j < 3; j++) {
139       crossProduct(axes[ct], axes[j], tri_edges[i]);
140       ct++;
141     }
142   }
143 }
144
145 /**
146  * Construction from a cube (axes aligned) and triangle
147  */
148 CubeTriangleIsect::CubeTriangleIsect(int64_t cube[2][3],
149                                      int64_t tri[3][3],
150                                      int64_t /*error*/,
151                                      int triind)
152 {
153   int i;
154   inherit = new TriangleProjection;
155   inherit->index = triind;
156
157   int64_t axes[NUM_AXES][3];
158   create_projection_axes(axes, tri);
159
160   /* Normalize face normal and store */
161   double dedge1[] = {(double)tri[1][0] - (double)tri[0][0],
162                      (double)tri[1][1] - (double)tri[0][1],
163                      (double)tri[1][2] - (double)tri[0][2]};
164   double dedge2[] = {(double)tri[2][0] - (double)tri[1][0],
165                      (double)tri[2][1] - (double)tri[1][1],
166                      (double)tri[2][2] - (double)tri[1][2]};
167   crossProduct(inherit->norm, dedge1, dedge2);
168   normalize(inherit->norm);
169
170   int64_t cubeedge[3][3];
171   for (i = 0; i < 3; i++) {
172     for (int j = 0; j < 3; j++) {
173       cubeedge[i][j] = 0;
174     }
175     cubeedge[i][i] = cube[1][i] - cube[0][i];
176   }
177
178   /* Project the cube on to each axis */
179   for (int axis = 0; axis < NUM_AXES; axis++) {
180     CubeProjection &cube_proj = cubeProj[axis];
181
182     /* Origin */
183     cube_proj.origin = dotProduct(axes[axis], cube[0]);
184
185     /* 3 direction vectors */
186     for (i = 0; i < 3; i++)
187       cube_proj.edges[i] = dotProduct(axes[axis], cubeedge[i]);
188
189     /* Offsets of 2 ends of cube projection */
190     int64_t max = 0;
191     int64_t min = 0;
192     for (i = 1; i < 8; i++) {
193       int64_t proj = (vertmap[i][0] * cube_proj.edges[0] + vertmap[i][1] * cube_proj.edges[1] +
194                       vertmap[i][2] * cube_proj.edges[2]);
195       if (proj > max) {
196         max = proj;
197       }
198       if (proj < min) {
199         min = proj;
200       }
201     }
202     cube_proj.min = min;
203     cube_proj.max = max;
204   }
205
206   /* Project the triangle on to each axis */
207   for (int axis = 0; axis < NUM_AXES; axis++) {
208     const int64_t vts[3] = {dotProduct(axes[axis], tri[0]),
209                             dotProduct(axes[axis], tri[1]),
210                             dotProduct(axes[axis], tri[2])};
211
212     // Triangle
213     inherit->tri_proj[axis][0] = vts[0];
214     inherit->tri_proj[axis][1] = vts[0];
215     for (i = 1; i < 3; i++) {
216       if (vts[i] < inherit->tri_proj[axis][0])
217         inherit->tri_proj[axis][0] = vts[i];
218
219       if (vts[i] > inherit->tri_proj[axis][1])
220         inherit->tri_proj[axis][1] = vts[i];
221     }
222   }
223 }
224
225 /**
226  * Construction
227  * from a parent CubeTriangleIsect object and the index of the children
228  */
229 CubeTriangleIsect::CubeTriangleIsect(CubeTriangleIsect *parent)
230 {
231   // Copy inheritable projections
232   this->inherit = parent->inherit;
233
234   // Shrink cube projections
235   for (int i = 0; i < NUM_AXES; i++) {
236     cubeProj[i].origin = parent->cubeProj[i].origin;
237
238     for (int j = 0; j < 3; j++)
239       cubeProj[i].edges[j] = parent->cubeProj[i].edges[j] >> 1;
240
241     cubeProj[i].min = parent->cubeProj[i].min >> 1;
242     cubeProj[i].max = parent->cubeProj[i].max >> 1;
243   }
244 }
245
246 unsigned char CubeTriangleIsect::getBoxMask()
247 {
248   int i, j, k;
249   int bmask[3][2] = {{0, 0}, {0, 0}, {0, 0}};
250   unsigned char boxmask = 0;
251   int64_t child_len = cubeProj[0].edges[0] >> 1;
252
253   for (i = 0; i < 3; i++) {
254     int64_t mid = cubeProj[i].origin + child_len;
255
256     // Check bounding box
257     if (mid >= inherit->tri_proj[i][0]) {
258       bmask[i][0] = 1;
259     }
260     if (mid < inherit->tri_proj[i][1]) {
261       bmask[i][1] = 1;
262     }
263   }
264
265   // Fill in masks
266   int ct = 0;
267   for (i = 0; i < 2; i++) {
268     for (j = 0; j < 2; j++) {
269       for (k = 0; k < 2; k++) {
270         boxmask |= ((bmask[0][i] & bmask[1][j] & bmask[2][k]) << ct);
271         ct++;
272       }
273     }
274   }
275
276   // Return bounding box masks
277   return boxmask;
278 }
279
280 /**
281  * Shifting a cube to a new origin
282  */
283 void CubeTriangleIsect::shift(int off[3])
284 {
285   for (int i = 0; i < NUM_AXES; i++) {
286     cubeProj[i].origin += (off[0] * cubeProj[i].edges[0] + off[1] * cubeProj[i].edges[1] +
287                            off[2] * cubeProj[i].edges[2]);
288   }
289 }
290
291 /**
292  * Method to test intersection of the triangle and the cube
293  */
294 int CubeTriangleIsect::isIntersecting() const
295 {
296   for (int i = 0; i < NUM_AXES; i++) {
297     /*
298       int64_t proj0 = cubeProj[i][0] +
299       vertmap[inherit->cubeEnds[i][0]][0] * cubeProj[i][1] +
300       vertmap[inherit->cubeEnds[i][0]][1] * cubeProj[i][2] +
301       vertmap[inherit->cubeEnds[i][0]][2] * cubeProj[i][3] ;
302       int64_t proj1 = cubeProj[i][0] +
303       vertmap[inherit->cubeEnds[i][1]][0] * cubeProj[i][1] +
304       vertmap[inherit->cubeEnds[i][1]][1] * cubeProj[i][2] +
305       vertmap[inherit->cubeEnds[i][1]][2] * cubeProj[i][3] ;
306     */
307
308     int64_t proj0 = cubeProj[i].origin + cubeProj[i].min;
309     int64_t proj1 = cubeProj[i].origin + cubeProj[i].max;
310
311     if (proj0 > inherit->tri_proj[i][1] || proj1 < inherit->tri_proj[i][0]) {
312       return 0;
313     }
314   }
315
316   return 1;
317 }
318
319 int CubeTriangleIsect::isIntersectingPrimary(int edgeInd) const
320 {
321   for (int i = 0; i < NUM_AXES; i++) {
322
323     int64_t proj0 = cubeProj[i].origin;
324     int64_t proj1 = cubeProj[i].origin + cubeProj[i].edges[edgeInd];
325
326     if (proj0 < proj1) {
327       if (proj0 > inherit->tri_proj[i][1] || proj1 < inherit->tri_proj[i][0]) {
328         return 0;
329       }
330     }
331     else {
332       if (proj1 > inherit->tri_proj[i][1] || proj0 < inherit->tri_proj[i][0]) {
333         return 0;
334       }
335     }
336   }
337
338   // printf( "Intersecting: %d %d\n", edgemap[edgeInd][0], edgemap[edgeInd][1] )  ;
339   return 1;
340 }
341
342 float CubeTriangleIsect::getIntersectionPrimary(int edgeInd) const
343 {
344   int i = 3;
345
346   int64_t proj0 = cubeProj[i].origin;
347   int64_t proj1 = cubeProj[i].origin + cubeProj[i].edges[edgeInd];
348   int64_t proj2 = inherit->tri_proj[i][1];
349   int64_t d = proj1 - proj0;
350   double alpha;
351
352   if (d == 0)
353     alpha = 0.5;
354   else {
355     alpha = (double)((proj2 - proj0)) / (double)d;
356
357     if (alpha < 0 || alpha > 1)
358       alpha = 0.5;
359   }
360
361   return (float)alpha;
362 }