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