Fix T39608: Blender 2.70 crashes when performing union
[blender.git] / extern / carve / carve-util.cc
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  * The Original Code is Copyright (C) 2014 Blender Foundation.
19  * All rights reserved.
20  *
21  * Contributor(s): Blender Foundation,
22  *                 Sergey Sharybin
23  *
24  * ***** END GPL LICENSE BLOCK *****
25  */
26
27 #include "carve-util.h"
28 #include "carve-capi.h"
29
30 #include <cfloat>
31 #include <carve/csg.hpp>
32 #include <carve/csg_triangulator.hpp>
33 #include <carve/rtree.hpp>
34
35 using carve::csg::Intersections;
36 using carve::geom::aabb;
37 using carve::geom::RTreeNode;
38 using carve::geom3d::Vector;
39 using carve::math::Matrix3;
40 using carve::mesh::Face;
41 using carve::mesh::MeshSet;
42 using carve::triangulate::tri_idx;
43 using carve::triangulate::triangulate;
44
45 typedef std::map< MeshSet<3>::mesh_t*, RTreeNode<3, Face<3> *> * > RTreeCache;
46 typedef std::map< MeshSet<3>::mesh_t*, bool > IntersectCache;
47
48 namespace {
49
50 // Functions adopted from BLI_math.h to use Carve Vector and Matrix.
51
52 void axis_angle_normalized_to_mat3(const Vector &normal,
53                                    const double angle,
54                                    Matrix3 *matrix)
55 {
56         double nsi[3], co, si, ico;
57
58         /* now convert this to a 3x3 matrix */
59         co = cos(angle);
60         si = sin(angle);
61
62         ico = (1.0 - co);
63         nsi[0] = normal[0] * si;
64         nsi[1] = normal[1] * si;
65         nsi[2] = normal[2] * si;
66
67         matrix->m[0][0] = ((normal[0] * normal[0]) * ico) + co;
68         matrix->m[0][1] = ((normal[0] * normal[1]) * ico) + nsi[2];
69         matrix->m[0][2] = ((normal[0] * normal[2]) * ico) - nsi[1];
70         matrix->m[1][0] = ((normal[0] * normal[1]) * ico) - nsi[2];
71         matrix->m[1][1] = ((normal[1] * normal[1]) * ico) + co;
72         matrix->m[1][2] = ((normal[1] * normal[2]) * ico) + nsi[0];
73         matrix->m[2][0] = ((normal[0] * normal[2]) * ico) + nsi[1];
74         matrix->m[2][1] = ((normal[1] * normal[2]) * ico) - nsi[0];
75         matrix->m[2][2] = ((normal[2] * normal[2]) * ico) + co;
76 }
77
78 void axis_angle_to_mat3(const Vector &axis,
79                         const double angle,
80                         Matrix3 *matrix)
81 {
82         if (axis.length2() < FLT_EPSILON) {
83                 *matrix = Matrix3();
84                 return;
85         }
86
87         Vector nor = axis;
88         nor.normalize();
89
90         axis_angle_normalized_to_mat3(nor, angle, matrix);
91 }
92
93 inline double saacos(double fac)
94 {
95         if      (fac <= -1.0) return M_PI;
96         else if (fac >=  1.0) return 0.0;
97         else                  return acos(fac);
98 }
99
100 bool axis_dominant_v3_to_m3(const Vector &normal,
101                             Matrix3 *matrix)
102 {
103         Vector up;
104         Vector axis;
105         double angle;
106
107         up.x = 0.0;
108         up.y = 0.0;
109         up.z = 1.0;
110
111         axis = carve::geom::cross(normal, up);
112         angle = saacos(carve::geom::dot(normal, up));
113
114         if (angle >= FLT_EPSILON) {
115                 if (axis.length2() < FLT_EPSILON) {
116                         axis[0] = 0.0;
117                         axis[1] = 1.0;
118                         axis[2] = 0.0;
119                 }
120
121                 axis_angle_to_mat3(axis, angle, matrix);
122                 return true;
123         }
124         else {
125                 *matrix = Matrix3();
126                 return false;
127         }
128 }
129
130 void meshset_minmax(const MeshSet<3> *mesh,
131                     Vector *min,
132                     Vector *max)
133 {
134         for (size_t i = 0; i < mesh->vertex_storage.size(); ++i) {
135                 min->x = std::min(min->x, mesh->vertex_storage[i].v.x);
136                 min->y = std::min(min->y, mesh->vertex_storage[i].v.y);
137                 min->z = std::min(min->z, mesh->vertex_storage[i].v.z);
138                 max->x = std::max(max->x, mesh->vertex_storage[i].v.x);
139                 max->y = std::max(max->y, mesh->vertex_storage[i].v.y);
140                 max->z = std::max(max->z, mesh->vertex_storage[i].v.z);
141         }
142 }
143
144 }  // namespace
145
146 void carve_getRescaleMinMax(const MeshSet<3> *left,
147                             const MeshSet<3> *right,
148                             Vector *min,
149                             Vector *max)
150 {
151         min->x = max->x = left->vertex_storage[0].v.x;
152         min->y = max->y = left->vertex_storage[0].v.y;
153         min->z = max->z = left->vertex_storage[0].v.z;
154
155         meshset_minmax(left, min, max);
156         meshset_minmax(right, min, max);
157
158         // Make sure we don't scale object with zero scale.
159         if (std::abs(min->x - max->x) < carve::EPSILON) {
160                 min->x = -1.0;
161                 max->x = 1.0;
162         }
163         if (std::abs(min->y - max->y) < carve::EPSILON) {
164                 min->y = -1.0;
165                 max->y = 1.0;
166         }
167         if (std::abs(min->z - max->z) < carve::EPSILON) {
168                 min->z = -1.0;
169                 max->z = 1.0;
170         }
171 }
172
173 namespace {
174
175 void copyMeshes(const std::vector<MeshSet<3>::mesh_t*> &meshes,
176                 std::vector<MeshSet<3>::mesh_t*> *new_meshes)
177 {
178         std::vector<MeshSet<3>::mesh_t*>::const_iterator it = meshes.begin();
179         new_meshes->reserve(meshes.size());
180         for (; it != meshes.end(); it++) {
181                 MeshSet<3>::mesh_t *mesh = *it;
182                 MeshSet<3>::mesh_t *new_mesh = new MeshSet<3>::mesh_t(mesh->faces);
183
184                 new_meshes->push_back(new_mesh);
185         }
186 }
187
188 MeshSet<3> *meshSetFromMeshes(const std::vector<MeshSet<3>::mesh_t*> &meshes)
189 {
190         std::vector<MeshSet<3>::mesh_t*> new_meshes;
191
192         copyMeshes(meshes, &new_meshes);
193
194         return new MeshSet<3>(new_meshes);
195 }
196
197 MeshSet<3> *meshSetFromTwoMeshes(const std::vector<MeshSet<3>::mesh_t*> &left_meshes,
198                                  const std::vector<MeshSet<3>::mesh_t*> &right_meshes)
199 {
200         std::vector<MeshSet<3>::mesh_t*> new_meshes;
201
202         copyMeshes(left_meshes, &new_meshes);
203         copyMeshes(right_meshes, &new_meshes);
204
205         return new MeshSet<3>(new_meshes);
206 }
207
208 bool checkEdgeFaceIntersections_do(Intersections &intersections,
209                                    MeshSet<3>::face_t *face_a,
210                                    MeshSet<3>::edge_t *edge_b)
211 {
212         if (intersections.intersects(edge_b, face_a))
213                 return true;
214
215         carve::mesh::MeshSet<3>::vertex_t::vector_t _p;
216         if (face_a->simpleLineSegmentIntersection(carve::geom3d::LineSegment(edge_b->v1()->v, edge_b->v2()->v), _p))
217                 return true;
218
219         return false;
220 }
221
222 bool checkEdgeFaceIntersections(Intersections &intersections,
223                                 MeshSet<3>::face_t *face_a,
224                                 MeshSet<3>::face_t *face_b)
225 {
226         MeshSet<3>::edge_t *edge_b;
227
228         edge_b = face_b->edge;
229         do {
230                 if (checkEdgeFaceIntersections_do(intersections, face_a, edge_b))
231                         return true;
232                 edge_b = edge_b->next;
233         } while (edge_b != face_b->edge);
234
235         return false;
236 }
237
238 inline bool facesAreCoplanar(const MeshSet<3>::face_t *a, const MeshSet<3>::face_t *b)
239 {
240         carve::geom3d::Ray temp;
241         // XXX: Find a better definition. This may be a source of problems
242         // if floating point inaccuracies cause an incorrect answer.
243         return !carve::geom3d::planeIntersection(a->plane, b->plane, temp);
244 }
245
246 bool checkMeshSetInterseciton_do(Intersections &intersections,
247                                  const RTreeNode<3, Face<3> *> *a_node,
248                                  const RTreeNode<3, Face<3> *> *b_node,
249                                  bool descend_a = true)
250 {
251         if (!a_node->bbox.intersects(b_node->bbox)) {
252                 return false;
253         }
254
255         if (a_node->child && (descend_a || !b_node->child)) {
256                 for (RTreeNode<3, Face<3> *> *node = a_node->child; node; node = node->sibling) {
257                         if (checkMeshSetInterseciton_do(intersections, node, b_node, false)) {
258                                 return true;
259                         }
260                 }
261         }
262         else if (b_node->child) {
263                 for (RTreeNode<3, Face<3> *> *node = b_node->child; node; node = node->sibling) {
264                         if (checkMeshSetInterseciton_do(intersections, a_node, node, true)) {
265                                 return true;
266                         }
267                 }
268         }
269         else {
270                 for (size_t i = 0; i < a_node->data.size(); ++i) {
271                         MeshSet<3>::face_t *fa = a_node->data[i];
272                         aabb<3> aabb_a = fa->getAABB();
273                         if (aabb_a.maxAxisSeparation(b_node->bbox) > carve::EPSILON) {
274                                 continue;
275                         }
276
277                         for (size_t j = 0; j < b_node->data.size(); ++j) {
278                                 MeshSet<3>::face_t *fb = b_node->data[j];
279                                 aabb<3> aabb_b = fb->getAABB();
280                                 if (aabb_b.maxAxisSeparation(aabb_a) > carve::EPSILON) {
281                                         continue;
282                                 }
283
284                                 std::pair<double, double> a_ra = fa->rangeInDirection(fa->plane.N, fa->edge->vert->v);
285                                 std::pair<double, double> b_ra = fb->rangeInDirection(fa->plane.N, fa->edge->vert->v);
286                                 if (carve::rangeSeparation(a_ra, b_ra) > carve::EPSILON) {
287                                         continue;
288                                 }
289
290                                 std::pair<double, double> a_rb = fa->rangeInDirection(fb->plane.N, fb->edge->vert->v);
291                                 std::pair<double, double> b_rb = fb->rangeInDirection(fb->plane.N, fb->edge->vert->v);
292                                 if (carve::rangeSeparation(a_rb, b_rb) > carve::EPSILON) {
293                                         continue;
294                                 }
295
296                                 if (!facesAreCoplanar(fa, fb)) {
297                                         if (checkEdgeFaceIntersections(intersections, fa, fb)) {
298                                                 return true;
299                                         }
300                                 }
301                         }
302                 }
303         }
304
305         return false;
306 }
307
308 bool checkMeshSetInterseciton(RTreeNode<3, Face<3> *> *rtree_a, RTreeNode<3, Face<3> *> *rtree_b)
309 {
310         Intersections intersections;
311         return checkMeshSetInterseciton_do(intersections, rtree_a, rtree_b);
312 }
313
314 void getIntersectedOperandMeshes(std::vector<MeshSet<3>::mesh_t*> *meshes,
315                                  const MeshSet<3>::aabb_t &otherAABB,
316                                  std::vector<MeshSet<3>::mesh_t*> *operandMeshes,
317                                  RTreeCache *rtree_cache,
318                                  IntersectCache *intersect_cache)
319 {
320         std::vector<MeshSet<3>::mesh_t*>::iterator it = meshes->begin();
321         std::vector< RTreeNode<3, Face<3> *> *> meshRTree;
322
323         while (it != meshes->end()) {
324                 MeshSet<3>::mesh_t *mesh = *it;
325                 bool isAdded = false;
326
327                 RTreeNode<3, Face<3> *> *rtree;
328                 bool intersects;
329
330                 RTreeCache::iterator rtree_found = rtree_cache->find(mesh);
331                 if (rtree_found != rtree_cache->end()) {
332                         rtree = rtree_found->second;
333                 }
334                 else {
335                         rtree = RTreeNode<3, Face<3> *>::construct_STR(mesh->faces.begin(), mesh->faces.end(), 4, 4);
336                         (*rtree_cache)[mesh] = rtree;
337                 }
338
339                 IntersectCache::iterator intersect_found = intersect_cache->find(mesh);
340                 if (intersect_found != intersect_cache->end()) {
341                         intersects = intersect_found->second;
342                 }
343                 else {
344                         intersects = rtree->bbox.intersects(otherAABB);
345                         (*intersect_cache)[mesh] = intersects;
346                 }
347
348                 if (intersects) {
349                         bool isIntersect = false;
350
351                         std::vector<MeshSet<3>::mesh_t*>::iterator operand_it = operandMeshes->begin();
352                         std::vector<RTreeNode<3, Face<3> *> *>::iterator tree_it = meshRTree.begin();
353                         for (; operand_it!=operandMeshes->end(); operand_it++, tree_it++) {
354                                 RTreeNode<3, Face<3> *> *operandRTree = *tree_it;
355
356                                 if (checkMeshSetInterseciton(rtree, operandRTree)) {
357                                         isIntersect = true;
358                                         break;
359                                 }
360                         }
361
362                         if (!isIntersect) {
363                                 operandMeshes->push_back(mesh);
364                                 meshRTree.push_back(rtree);
365
366                                 it = meshes->erase(it);
367                                 isAdded = true;
368                         }
369                 }
370
371                 if (!isAdded) {
372                         //delete rtree;
373                         it++;
374                 }
375         }
376
377         std::vector<RTreeNode<3, Face<3> *> *>::iterator tree_it = meshRTree.begin();
378         for (; tree_it != meshRTree.end(); tree_it++) {
379                 //delete *tree_it;
380         }
381 }
382
383 MeshSet<3> *getIntersectedOperand(std::vector<MeshSet<3>::mesh_t*> *meshes,
384                                   const MeshSet<3>::aabb_t &otherAABB,
385                                   RTreeCache *rtree_cache,
386                                   IntersectCache *intersect_cache)
387 {
388         std::vector<MeshSet<3>::mesh_t*> operandMeshes;
389         getIntersectedOperandMeshes(meshes, otherAABB, &operandMeshes, rtree_cache, intersect_cache);
390
391         if (operandMeshes.size() == 0)
392                 return NULL;
393
394         return meshSetFromMeshes(operandMeshes);
395 }
396
397 MeshSet<3> *unionIntersectingMeshes(carve::csg::CSG *csg,
398                                     MeshSet<3> *poly,
399                                     const MeshSet<3>::aabb_t &otherAABB)
400 {
401         if (poly->meshes.size() <= 1) {
402                 return poly;
403         }
404
405         std::vector<MeshSet<3>::mesh_t*> orig_meshes =
406                         std::vector<MeshSet<3>::mesh_t*>(poly->meshes.begin(), poly->meshes.end());
407
408         RTreeCache rtree_cache;
409         IntersectCache intersect_cache;
410
411         MeshSet<3> *left = getIntersectedOperand(&orig_meshes,
412                                                  otherAABB,
413                                                  &rtree_cache,
414                                                  &intersect_cache);
415
416         if (!left) {
417                 // No maniforlds which intersects another object at all.
418                 return poly;
419         }
420
421         while (orig_meshes.size()) {
422                 MeshSet<3> *right = getIntersectedOperand(&orig_meshes,
423                                                           otherAABB,
424                                                           &rtree_cache,
425                                                           &intersect_cache);
426
427                 if (!right) {
428                         // No more intersecting manifolds which intersects other object
429                         break;
430                 }
431
432                 try {
433                         if (left->meshes.size()==0) {
434                                 delete left;
435
436                                 left = right;
437                         }
438                         else {
439                                 MeshSet<3> *result = csg->compute(left, right,
440                                                                   carve::csg::CSG::UNION,
441                                                                   NULL, carve::csg::CSG::CLASSIFY_EDGE);
442
443                                 delete left;
444                                 delete right;
445
446                                 left = result;
447                         }
448                 }
449                 catch (carve::exception e) {
450                         std::cerr << "CSG failed, exception " << e.str() << std::endl;
451
452                         MeshSet<3> *result = meshSetFromTwoMeshes(left->meshes, right->meshes);
453
454                         delete left;
455                         delete right;
456
457                         left = result;
458                 }
459                 catch (...) {
460                         delete left;
461                         delete right;
462
463                         throw "Unknown error in Carve library";
464                 }
465         }
466
467         for (RTreeCache::iterator it = rtree_cache.begin();
468              it != rtree_cache.end();
469              it++)
470         {
471                 delete it->second;
472         }
473
474         // Append all meshes which doesn't have intersection with another operand as-is.
475         if (orig_meshes.size()) {
476                 MeshSet<3> *result = meshSetFromTwoMeshes(left->meshes, orig_meshes);
477
478                 delete left;
479                 left = result;
480         }
481
482         return left;
483 }
484
485 }  // namespace
486
487 // TODO(sergey): This function is to be totally re-implemented to make it
488 // more clear what's going on and hopefully optimize it as well.
489 bool carve_unionIntersections(carve::csg::CSG *csg,
490                               MeshSet<3> **left_r,
491                               MeshSet<3> **right_r)
492 {
493         MeshSet<3> *left = *left_r, *right = *right_r;
494         bool changed = false;
495
496         if (left->meshes.size() == 1 && right->meshes.size() == 0) {
497                 return false;
498         }
499
500         MeshSet<3>::aabb_t leftAABB = left->getAABB();
501         MeshSet<3>::aabb_t rightAABB = right->getAABB();;
502
503         left = unionIntersectingMeshes(csg, left, rightAABB);
504         right = unionIntersectingMeshes(csg, right, leftAABB);
505
506         if (left != *left_r) {
507                 changed = true;
508                 delete *left_r;
509         }
510
511         if (right != *right_r) {
512                 changed = true;
513                 delete *right_r;
514         }
515
516         *left_r = left;
517         *right_r = right;
518
519         return changed;
520 }
521
522 static inline void add_newell_cross_v3_v3v3(const Vector &v_prev,
523                                             const Vector &v_curr,
524                                             Vector *n)
525 {
526         (*n)[0] += (v_prev[1] - v_curr[1]) * (v_prev[2] + v_curr[2]);
527         (*n)[1] += (v_prev[2] - v_curr[2]) * (v_prev[0] + v_curr[0]);
528         (*n)[2] += (v_prev[0] - v_curr[0]) * (v_prev[1] + v_curr[1]);
529 }
530
531 // Axis matrix is being set for non-flat ngons only.
532 bool carve_checkPolyPlanarAndGetNormal(const std::vector<Vector> &vertices,
533                                        const int verts_per_poly,
534                                        const int *verts_of_poly,
535                                        Matrix3 *axis_matrix_r)
536 {
537         if (verts_per_poly == 3) {
538                 // Triangles are always planar.
539                 return true;
540         }
541         else if (verts_per_poly == 4) {
542                 // Presumably faster than using generig n-gon check for quads.
543
544                 const Vector &v1 = vertices[verts_of_poly[0]],
545                              &v2 = vertices[verts_of_poly[1]],
546                              &v3 = vertices[verts_of_poly[2]],
547                              &v4 = vertices[verts_of_poly[3]];
548
549                 Vector vec1, vec2, vec3, cross;
550
551                 vec1 = v2 - v1;
552                 vec2 = v4 - v1;
553                 vec3 = v3 - v1;
554
555                 cross = carve::geom::cross(vec1, vec2);
556
557                 double production = carve::geom::dot(cross, vec3);
558
559                 // TODO(sergey): Check on whether we could have length-independent
560                 // magnitude here.
561                 double magnitude = 1e-3 * cross.length2();
562
563                 return fabs(production) < magnitude;
564         }
565         else {
566                 const Vector *vert_prev = &vertices[verts_of_poly[verts_per_poly - 1]];
567                 const Vector *vert_curr = &vertices[verts_of_poly[0]];
568
569                 Vector normal = carve::geom::VECTOR(0.0, 0.0, 0.0);
570                 for (int i = 0; i < verts_per_poly; i++) {
571                         add_newell_cross_v3_v3v3(*vert_prev, *vert_curr, &normal);
572                         vert_prev = vert_curr;
573                         vert_curr = &vertices[verts_of_poly[(i + 1) % verts_per_poly]];
574                 }
575
576                 if (normal.length2() < FLT_EPSILON) {
577                         // Degenerated face, couldn't triangulate properly anyway.
578                         return true;
579                 }
580                 else {
581                         double magnitude = normal.length2();
582
583                         normal.normalize();
584                         axis_dominant_v3_to_m3(normal, axis_matrix_r);
585
586                         Vector first_projected = *axis_matrix_r * vertices[verts_of_poly[0]];
587                         double min_z = first_projected[2], max_z = first_projected[2];
588
589                         for (int i = 1; i < verts_per_poly; i++) {
590                                 const Vector &vertex = vertices[verts_of_poly[i]];
591                                 Vector projected = *axis_matrix_r * vertex;
592                                 if (projected[2] < min_z) {
593                                         min_z = projected[2];
594                                 }
595                                 if (projected[2] > max_z) {
596                                         max_z = projected[2];
597                                 }
598                         }
599
600                         if (std::abs(min_z - max_z) > FLT_EPSILON * magnitude) {
601                                 return false;
602                         }
603                 }
604
605                 return true;
606         }
607
608         return false;
609 }
610
611 namespace {
612
613 int triangulateNGon_carveTriangulator(const std::vector<Vector> &vertices,
614                                       const int verts_per_poly,
615                                       const int *verts_of_poly,
616                                       const Matrix3 &axis_matrix,
617                                       std::vector<tri_idx> *triangles)
618 {
619         // Project vertices to 2D plane.
620         Vector projected;
621         std::vector<carve::geom::vector<2> > poly_2d;
622         poly_2d.reserve(verts_per_poly);
623         for (int i = 0; i < verts_per_poly; ++i) {
624                 projected = axis_matrix * vertices[verts_of_poly[i]];
625                 poly_2d.push_back(carve::geom::VECTOR(projected[0], projected[1]));
626         }
627
628         carve::triangulate::triangulate(poly_2d, *triangles);
629         carve::triangulate::improve(poly_2d, *triangles);
630
631         return triangles->size();
632 }
633
634 int triangulateNGon_importerTriangulator(struct ImportMeshData *import_data,
635                                          CarveMeshImporter *mesh_importer,
636                                          const std::vector<Vector> &vertices,
637                                          const int verts_per_poly,
638                                          const int *verts_of_poly,
639                                          const Matrix3 &axis_matrix,
640                                          std::vector<tri_idx> *triangles)
641 {
642         typedef float Vector2D[2];
643         typedef unsigned int Triangle[3];
644
645         // Project vertices to 2D plane.
646         Vector2D *poly_2d = new Vector2D[verts_per_poly];
647         Vector projected;
648         for (int i = 0; i < verts_per_poly; ++i) {
649                 projected = axis_matrix * vertices[verts_of_poly[i]];
650                 poly_2d[i][0] = projected[0];
651                 poly_2d[i][1] = projected[1];
652         }
653
654         Triangle *api_triangles = new Triangle[verts_per_poly - 2];
655         int num_triangles =
656                 mesh_importer->triangulate2DPoly(import_data,
657                                                  poly_2d,
658                                                  verts_per_poly,
659                                                  api_triangles);
660
661         triangles->reserve(num_triangles);
662         for (int i = 0; i < num_triangles; ++i) {
663                 triangles->push_back(tri_idx(api_triangles[i][0],
664                                                  api_triangles[i][1],
665                                                  api_triangles[i][2]));
666         }
667
668         delete [] poly_2d;
669         delete [] api_triangles;
670
671         return num_triangles;
672 }
673
674 template <typename T>
675 void sortThreeNumbers(T &a, T &b, T &c)
676 {
677         if (a > b)
678                 std::swap(a, b);
679         if (b > c)
680                 std::swap(b, c);
681         if (a > b)
682                 std::swap(a, b);
683 }
684
685 bool pushTriangle(int v1, int v2, int v3,
686                   std::vector<int> *face_indices,
687                   TrianglesStorage *triangles_storage)
688 {
689
690         tri_idx triangle(v1, v2, v3);
691         sortThreeNumbers(triangle.a, triangle.b, triangle.c);
692
693         assert(triangle.a < triangle.b);
694         assert(triangle.b < triangle.c);
695
696         if (triangles_storage->find(triangle) == triangles_storage->end()) {
697                 face_indices->push_back(3);
698                 face_indices->push_back(v1);
699                 face_indices->push_back(v2);
700                 face_indices->push_back(v3);
701
702                 triangles_storage->insert(triangle);
703                 return true;
704         }
705         else {
706                 return false;
707         }
708 }
709
710 }  // namespace
711
712 int carve_triangulatePoly(struct ImportMeshData *import_data,
713                           CarveMeshImporter *mesh_importer,
714                           const std::vector<Vector> &vertices,
715                           const int verts_per_poly,
716                           const int *verts_of_poly,
717                           const Matrix3 &axis_matrix,
718                           std::vector<int> *face_indices,
719                           TrianglesStorage *triangles_storage)
720 {
721         int num_triangles = 0;
722
723         assert(verts_per_poly > 3);
724
725         if (verts_per_poly == 4) {
726                 // Quads we triangulate by 1-3 diagonal, it is an original behavior
727                 // of boolean modifier.
728                 //
729                 // TODO(sergey): Consider using shortest diagonal here. However
730                 // display code in Blende use static 1-3 split, so some experiments
731                 // are needed here.
732                 if (pushTriangle(verts_of_poly[0],
733                                  verts_of_poly[1],
734                                  verts_of_poly[2],
735                                  face_indices,
736                                  triangles_storage))
737                 {
738                         num_triangles++;
739                 }
740
741                 if (pushTriangle(verts_of_poly[0],
742                                  verts_of_poly[2],
743                                  verts_of_poly[3],
744                                  face_indices,
745                                  triangles_storage))
746                 {
747                         num_triangles++;
748                 }
749         }
750         else {
751                 std::vector<tri_idx> triangles;
752                 triangles.reserve(verts_per_poly - 2);
753
754                 // Make triangulator callback optional so we could do some tests
755                 // in the future.
756                 if (mesh_importer->triangulate2DPoly) {
757                         triangulateNGon_importerTriangulator(import_data,
758                                                              mesh_importer,
759                                                              vertices,
760                                                              verts_per_poly,
761                                                              verts_of_poly,
762                                                              axis_matrix,
763                                                              &triangles);
764                 }
765                 else {
766                         triangulateNGon_carveTriangulator(vertices,
767                                                           verts_per_poly,
768                                                           verts_of_poly,
769                                                           axis_matrix,
770                                                           &triangles);
771                 }
772
773                 for (int i = 0; i < triangles.size(); ++i) {
774                         int v1 = triangles[i].a,
775                             v2 = triangles[i].b,
776                             v3 = triangles[i].c;
777
778                         // Sanity check of the triangle.
779                         assert(v1 != v2);
780                         assert(v1 != v3);
781                         assert(v2 != v3);
782                         assert(v1 < verts_per_poly);
783                         assert(v2 < verts_per_poly);
784                         assert(v3 < verts_per_poly);
785
786                         if (pushTriangle(verts_of_poly[v1],
787                                          verts_of_poly[v2],
788                                          verts_of_poly[v3],
789                                          face_indices,
790                                          triangles_storage))
791                         {
792                                 num_triangles++;
793                         }
794                 }
795         }
796
797         return num_triangles;
798 }