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