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