svn merge ^/trunk/blender -r49082:49104
[blender.git] / intern / boolop / intern / BOP_CarveInterface.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  * The Original Code is Copyright (C) 2010 by the Blender Foundation.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): Ken Hughes,
24  *                 Sergey Sharybin.
25  *
26  * ***** END GPL LICENSE BLOCK *****
27  */
28
29 /** \file boolop/intern/BOP_CarveInterface.cpp
30  *  \ingroup boolopintern
31  */
32
33 #include "../extern/BOP_Interface.h"
34 #include "../../bsp/intern/BSP_CSGMesh_CFIterator.h"
35
36 #include <carve/csg_triangulator.hpp>
37 #include <carve/interpolator.hpp>
38 #include <carve/rescale.hpp>
39
40 #include <iostream>
41
42 using namespace carve::mesh;
43 using namespace carve::geom;
44 typedef unsigned int uint;
45
46 #define MAX(x,y) ((x)>(y)?(x):(y))
47 #define MIN(x,y) ((x)<(y)?(x):(y))
48
49 static bool isQuadPlanar(carve::geom3d::Vector &v1, carve::geom3d::Vector &v2,
50                          carve::geom3d::Vector &v3, carve::geom3d::Vector &v4)
51 {
52         carve::geom3d::Vector vec1, vec2, vec3, cross;
53
54         vec1 = v2 - v1;
55         vec2 = v4 - v1;
56         vec3 = v3 - v1;
57
58         cross = carve::geom::cross(vec1, vec2);
59
60         float production = carve::geom::dot(cross, vec3);
61         float magnitude = 1e-6 * cross.length();
62
63         return fabs(production) < magnitude;
64 }
65
66 static bool isFacePlanar(CSG_IFace &face, std::vector<carve::geom3d::Vector> &vertices)
67 {
68         if (face.vertex_number == 4) {
69                 return isQuadPlanar(vertices[face.vertex_index[0]], vertices[face.vertex_index[1]],
70                                     vertices[face.vertex_index[2]], vertices[face.vertex_index[3]]);
71         }
72
73         return true;
74 }
75
76 static void Carve_copyMeshes(std::vector<MeshSet<3>::mesh_t*> &meshes, std::vector<MeshSet<3>::mesh_t*> &new_meshes)
77 {
78         std::vector<MeshSet<3>::mesh_t*>::iterator it = meshes.begin();
79
80         for(; it!=meshes.end(); it++) {
81                 MeshSet<3>::mesh_t *mesh = *it;
82                 MeshSet<3>::mesh_t *new_mesh = new MeshSet<3>::mesh_t(mesh->faces);
83
84                 new_meshes.push_back(new_mesh);
85         }
86 }
87
88 static MeshSet<3> *Carve_meshSetFromMeshes(std::vector<MeshSet<3>::mesh_t*> &meshes)
89 {
90         std::vector<MeshSet<3>::mesh_t*> new_meshes;
91
92         Carve_copyMeshes(meshes, new_meshes);
93
94         return new MeshSet<3>(new_meshes);
95 }
96
97 static MeshSet<3> *Carve_meshSetFromTwoMeshes(std::vector<MeshSet<3>::mesh_t*> &left_meshes,
98                                               std::vector<MeshSet<3>::mesh_t*> &right_meshes)
99 {
100         std::vector<MeshSet<3>::mesh_t*> new_meshes;
101
102         Carve_copyMeshes(left_meshes, new_meshes);
103         Carve_copyMeshes(right_meshes, new_meshes);
104
105         return new MeshSet<3>(new_meshes);
106 }
107
108 static bool Carve_checkEdgeFaceIntersections_do(carve::csg::Intersections &intersections,
109                                                 MeshSet<3>::face_t *face_a, MeshSet<3>::edge_t *edge_b)
110 {
111         if(intersections.intersects(edge_b, face_a))
112                 return true;
113
114         carve::mesh::MeshSet<3>::vertex_t::vector_t _p;
115         if(face_a->simpleLineSegmentIntersection(carve::geom3d::LineSegment(edge_b->v1()->v, edge_b->v2()->v), _p))
116                 return true;
117
118         return false;
119 }
120
121 static bool Carve_checkEdgeFaceIntersections(carve::csg::Intersections &intersections,
122                                              MeshSet<3>::face_t *face_a, MeshSet<3>::face_t *face_b)
123 {
124         MeshSet<3>::edge_t *edge_b;
125
126         edge_b = face_b->edge;
127         do {
128                 if(Carve_checkEdgeFaceIntersections_do(intersections, face_a, edge_b))
129                         return true;
130                 edge_b = edge_b->next;
131         } while (edge_b != face_b->edge);
132
133         return false;
134 }
135
136 static inline bool Carve_facesAreCoplanar(const MeshSet<3>::face_t *a, const MeshSet<3>::face_t *b)
137 {
138   carve::geom3d::Ray temp;
139   // XXX: Find a better definition. This may be a source of problems
140   // if floating point inaccuracies cause an incorrect answer.
141   return !carve::geom3d::planeIntersection(a->plane, b->plane, temp);
142 }
143
144 static bool Carve_checkMeshSetInterseciton_do(carve::csg::Intersections &intersections,
145                                               const RTreeNode<3, Face<3> *> *a_node,
146                                               const RTreeNode<3, Face<3> *> *b_node,
147                                               bool descend_a = true)
148 {
149         if(!a_node->bbox.intersects(b_node->bbox))
150                 return false;
151
152         if(a_node->child && (descend_a || !b_node->child)) {
153                 for(RTreeNode<3, Face<3> *> *node = a_node->child; node; node = node->sibling) {
154                         if(Carve_checkMeshSetInterseciton_do(intersections, node, b_node, false))
155                                 return true;
156                 }
157         }
158         else if(b_node->child) {
159                 for(RTreeNode<3, Face<3> *> *node = b_node->child; node; node = node->sibling) {
160                         if(Carve_checkMeshSetInterseciton_do(intersections, a_node, node, true))
161                                 return true;
162                 }
163         }
164         else {
165                 for(size_t i = 0; i < a_node->data.size(); ++i) {
166                         MeshSet<3>::face_t *fa = a_node->data[i];
167                         aabb<3> aabb_a = fa->getAABB();
168                         if(aabb_a.maxAxisSeparation(b_node->bbox) > carve::EPSILON) continue;
169
170                         for(size_t j = 0; j < b_node->data.size(); ++j) {
171                                 MeshSet<3>::face_t *fb = b_node->data[j];
172                                 aabb<3> aabb_b = fb->getAABB();
173                                 if(aabb_b.maxAxisSeparation(aabb_a) > carve::EPSILON) continue;
174
175                                 std::pair<double, double> a_ra = fa->rangeInDirection(fa->plane.N, fa->edge->vert->v);
176                                 std::pair<double, double> b_ra = fb->rangeInDirection(fa->plane.N, fa->edge->vert->v);
177                                 if(carve::rangeSeparation(a_ra, b_ra) > carve::EPSILON) continue;
178
179                                 std::pair<double, double> a_rb = fa->rangeInDirection(fb->plane.N, fb->edge->vert->v);
180                                 std::pair<double, double> b_rb = fb->rangeInDirection(fb->plane.N, fb->edge->vert->v);
181                                 if(carve::rangeSeparation(a_rb, b_rb) > carve::EPSILON) continue;
182
183                                 if(!Carve_facesAreCoplanar(fa, fb)) {
184                                         if(Carve_checkEdgeFaceIntersections(intersections, fa, fb)) {
185                                                 return true;
186                                         }
187                                 }
188                         }
189                 }
190         }
191
192         return false;
193 }
194
195 static bool Carve_checkMeshSetInterseciton(RTreeNode<3, Face<3> *> *rtree_a, RTreeNode<3, Face<3> *> *rtree_b)
196 {
197         carve::csg::Intersections intersections;
198
199         return Carve_checkMeshSetInterseciton_do(intersections, rtree_a, rtree_b);
200 }
201
202 static void Carve_getIntersectedOperandMeshes(std::vector<MeshSet<3>::mesh_t*> &meshes, MeshSet<3>::aabb_t &otherAABB,
203                                               std::vector<MeshSet<3>::mesh_t*> &operandMeshes)
204 {
205         std::vector<MeshSet<3>::mesh_t*>::iterator it = meshes.begin();
206         std::vector< RTreeNode<3, Face<3> *> *> meshRTree;
207
208         while(it != meshes.end()) {
209                 MeshSet<3>::mesh_t *mesh = *it;
210                 bool isAdded = false;
211
212                 RTreeNode<3, Face<3> *> *rtree = RTreeNode<3, Face<3> *>::construct_STR(mesh->faces.begin(), mesh->faces.end(), 4, 4);
213
214                 if (rtree->bbox.intersects(otherAABB)) {
215                         bool isIntersect = false;
216
217                         std::vector<MeshSet<3>::mesh_t*>::iterator operand_it = operandMeshes.begin();
218                         std::vector<RTreeNode<3, Face<3> *> *>::iterator tree_it = meshRTree.begin();
219                         for(; operand_it!=operandMeshes.end(); operand_it++, tree_it++) {
220                                 RTreeNode<3, Face<3> *> *operandRTree = *tree_it;
221
222                                 if(Carve_checkMeshSetInterseciton(rtree, operandRTree)) {
223                                         isIntersect = true;
224                                         break;
225                                 }
226                         }
227
228                         if(!isIntersect) {
229                                 operandMeshes.push_back(mesh);
230                                 meshRTree.push_back(rtree);
231
232                                 it = meshes.erase(it);
233                                 isAdded = true;
234                         }
235                 }
236
237                 if (!isAdded) {
238                         delete rtree;
239                         it++;
240                 }
241         }
242
243         std::vector<RTreeNode<3, Face<3> *> *>::iterator tree_it = meshRTree.begin();
244         for(; tree_it != meshRTree.end(); tree_it++) {
245                 delete *tree_it;
246         }
247 }
248
249 static MeshSet<3> *Carve_getIntersectedOperand(std::vector<MeshSet<3>::mesh_t*> &meshes, MeshSet<3>::aabb_t &otherAABB)
250 {
251         std::vector<MeshSet<3>::mesh_t*> operandMeshes;
252         Carve_getIntersectedOperandMeshes(meshes, otherAABB, operandMeshes);
253
254         if (operandMeshes.size() == 0)
255                 return NULL;
256
257         return Carve_meshSetFromMeshes(operandMeshes);
258 }
259
260 static MeshSet<3> *Carve_unionIntersectingMeshes(MeshSet<3> *poly,
261                                                  MeshSet<3>::aabb_t &otherAABB,
262                                                  carve::interpolate::FaceAttr<uint> &oface_num)
263 {
264         if(poly->meshes.size()<=1)
265                 return poly;
266
267         carve::csg::CSG csg;
268
269         oface_num.installHooks(csg);
270         csg.hooks.registerHook(new carve::csg::CarveTriangulator, carve::csg::CSG::Hooks::PROCESS_OUTPUT_FACE_BIT);
271
272         std::vector<MeshSet<3>::mesh_t*> orig_meshes =
273                         std::vector<MeshSet<3>::mesh_t*>(poly->meshes.begin(), poly->meshes.end());
274
275         MeshSet<3> *left = Carve_getIntersectedOperand(orig_meshes, otherAABB);
276
277         if (!left) {
278                 /* no maniforlds which intersects another object at all */
279                 return poly;
280         }
281
282         while(orig_meshes.size()) {
283                 MeshSet<3> *right = Carve_getIntersectedOperand(orig_meshes, otherAABB);
284
285                 if (!right) {
286                         /* no more intersecting manifolds which intersects other object */
287                         break;
288                 }
289
290                 try {
291                         if(left->meshes.size()==0) {
292                                 delete left;
293
294                                 left = right;
295                         }
296                         else {
297                                 MeshSet<3> *result = csg.compute(left, right, carve::csg::CSG::UNION, NULL, carve::csg::CSG::CLASSIFY_EDGE);
298
299                                 delete left;
300                                 delete right;
301
302                                 left = result;
303                         }
304                 }
305                 catch(carve::exception e) {
306                         std::cerr << "CSG failed, exception " << e.str() << std::endl;
307
308                         MeshSet<3> *result = Carve_meshSetFromTwoMeshes(left->meshes, right->meshes);
309
310                         delete left;
311                         delete right;
312
313                         left = result;
314                 }
315                 catch(...) {
316                         delete left;
317                         delete right;
318
319                         throw "Unknown error in Carve library";
320                 }
321         }
322
323         /* append all meshes which doesn't have intersection with another operand as-is */
324         if (orig_meshes.size()) {
325                 MeshSet<3> *result = Carve_meshSetFromTwoMeshes(left->meshes, orig_meshes);
326
327                 delete left;
328
329                 return result;
330         }
331
332         return left;
333 }
334
335 static void Carve_unionIntersections(MeshSet<3> **left_r, MeshSet<3> **right_r,
336                                      carve::interpolate::FaceAttr<uint> &oface_num)
337 {
338         MeshSet<3> *left, *right;
339
340         MeshSet<3>::aabb_t leftAABB = (*left_r)->getAABB();
341         MeshSet<3>::aabb_t rightAABB = (*right_r)->getAABB();
342
343         left = Carve_unionIntersectingMeshes(*left_r, rightAABB, oface_num);
344         right = Carve_unionIntersectingMeshes(*right_r, leftAABB, oface_num);
345
346         if(left != *left_r)
347                 delete *left_r;
348
349         if(right != *right_r)
350                 delete *right_r;
351
352         *left_r = left;
353         *right_r = right;
354 }
355
356 static MeshSet<3> *Carve_addMesh(CSG_FaceIteratorDescriptor &face_it,
357                                  CSG_VertexIteratorDescriptor &vertex_it,
358                                  carve::interpolate::FaceAttr<uint> &oface_num,
359                                  uint &num_origfaces)
360 {
361         CSG_IVertex vertex;
362         std::vector<carve::geom3d::Vector> vertices;
363
364         while (!vertex_it.Done(vertex_it.it)) {
365                 vertex_it.Fill(vertex_it.it,&vertex);
366                 vertices.push_back(VECTOR(vertex.position[0],
367                                           vertex.position[1],
368                                           vertex.position[2]));
369                 vertex_it.Step(vertex_it.it);
370         }
371
372         CSG_IFace face;
373         std::vector<int> f;
374         int numfaces = 0;
375
376         // now for the polygons.
377         // we may need to decalare some memory for user defined face properties.
378
379         std::vector<int> forig;
380         while (!face_it.Done(face_it.it)) {
381                 face_it.Fill(face_it.it,&face);
382
383                 if (isFacePlanar(face, vertices)) {
384                         f.push_back(face.vertex_number);
385                         f.push_back(face.vertex_index[0]);
386                         f.push_back(face.vertex_index[1]);
387                         f.push_back(face.vertex_index[2]);
388
389                         if (face.vertex_number == 4)
390                                 f.push_back(face.vertex_index[3]);
391
392                         forig.push_back(face.orig_face);
393                         ++numfaces;
394                         face_it.Step(face_it.it);
395                         ++num_origfaces;
396                 }
397                 else {
398                         f.push_back(3);
399                         f.push_back(face.vertex_index[0]);
400                         f.push_back(face.vertex_index[1]);
401                         f.push_back(face.vertex_index[2]);
402
403                         forig.push_back(face.orig_face);
404                         ++numfaces;
405
406                         if (face.vertex_number == 4) {
407                                 f.push_back(3);
408                                 f.push_back(face.vertex_index[0]);
409                                 f.push_back(face.vertex_index[2]);
410                                 f.push_back(face.vertex_index[3]);
411
412                                 forig.push_back(face.orig_face);
413                                 ++numfaces;
414                         }
415
416                         face_it.Step(face_it.it);
417                         ++num_origfaces;
418                 }
419         }
420
421         MeshSet<3> *poly = new MeshSet<3> (vertices, numfaces, f);
422
423         uint i;
424         MeshSet<3>::face_iter face_iter = poly->faceBegin();
425         for (i = 0; face_iter != poly->faceEnd(); ++face_iter, ++i) {
426                 MeshSet<3>::face_t *face = *face_iter;
427                 oface_num.setAttribute(face, forig[i]);
428         }
429
430         return poly;
431 }
432
433 static double triangleArea(carve::geom3d::Vector &v1, carve::geom3d::Vector &v2, carve::geom3d::Vector &v3)
434 {
435         carve::geom3d::Vector a = v2 - v1;
436         carve::geom3d::Vector b = v3 - v1;
437
438         return carve::geom::cross(a, b).length();
439 }
440
441 static bool checkValidQuad(std::vector<MeshSet<3>::vertex_t> &vertex_storage, uint quad[4])
442 {
443         carve::geom3d::Vector &v1 = vertex_storage[quad[0]].v;
444         carve::geom3d::Vector &v2 = vertex_storage[quad[1]].v;
445         carve::geom3d::Vector &v3 = vertex_storage[quad[2]].v;
446         carve::geom3d::Vector &v4 = vertex_storage[quad[3]].v;
447
448 #if 0
449         /* disabled for now to prevent initially non-planar be triangulated
450          * in theory this might cause some artifacts if intersections happens by non-planar
451          * non-concave quad, but in practice it's acceptable */
452         if (!isQuadPlanar(v1, v2, v3, v4)) {
453                 /* non-planar faces better not be merged because of possible differences in triangulation
454                  * of non-planar faces in opengl and renderer */
455                 return false;
456         }
457 #endif
458
459         carve::geom3d::Vector edges[4];
460         carve::geom3d::Vector normal;
461         bool normal_set = false;
462
463         edges[0] = v2 - v1;
464         edges[1] = v3 - v2;
465         edges[2] = v4 - v3;
466         edges[3] = v1 - v4;
467
468         for (int i = 0; i < 4; i++) {
469                 int n = i + 1;
470
471                 if (n == 4)
472                         n = 0;
473
474                 carve::geom3d::Vector current_normal = carve::geom::cross(edges[i], edges[n]);
475
476                 if (current_normal.length() > DBL_EPSILON) {
477                         if (!normal_set) {
478                                 normal = current_normal;
479                                 normal_set = true;
480                         }
481                         else if (carve::geom::dot(normal, current_normal) < 0) {
482                                 return false;
483                         }
484                 }
485         }
486
487         if (!normal_set) {
488                 /* normal wasn't set means face is degraded and better merge it in such way */
489                 return false;
490         }
491
492         double area = triangleArea(v1, v2, v3) + triangleArea(v1, v3, v4);
493         if (area <= DBL_EPSILON)
494                 return false;
495
496         return true;
497 }
498
499 // check whether two faces share an edge, and if so merge them
500 static uint quadMerge(std::map<MeshSet<3>::vertex_t*, uint> *vertexToIndex_map,
501                                           std::vector<MeshSet<3>::vertex_t> &vertex_storage,
502                       MeshSet<3>::face_t *f1, MeshSet<3>::face_t *f2,
503                       uint v, uint quad[4])
504 {
505         uint current, n1, p1, n2, p2;
506         uint v1[3];
507         uint v2[3];
508
509         // get the vertex indices for each face
510         v1[0] = vertexToIndex_map->find(f1->edge->vert)->second;
511         v1[1] = vertexToIndex_map->find(f1->edge->next->vert)->second;
512         v1[2] = vertexToIndex_map->find(f1->edge->next->next->vert)->second;
513
514         v2[0] = vertexToIndex_map->find(f2->edge->vert)->second;
515         v2[1] = vertexToIndex_map->find(f2->edge->next->vert)->second;
516         v2[2] = vertexToIndex_map->find(f2->edge->next->next->vert)->second;
517
518         // locate the current vertex we're examining, and find the next and
519         // previous vertices based on the face windings
520         if (v1[0] == v) {current = 0; p1 = 2; n1 = 1;}
521         else if (v1[1] == v) {current = 1; p1 = 0; n1 = 2;}
522         else {current = 2; p1 = 1; n1 = 0;}
523
524         if (v2[0] == v) {p2 = 2; n2 = 1;}
525         else if (v2[1] == v) {p2 = 0; n2 = 2;}
526         else {p2 = 1; n2 = 0;}
527
528         // if we find a match, place indices into quad in proper order and return
529         // success code
530         if (v1[p1] == v2[n2]) {
531                 quad[0] = v1[current];
532                 quad[1] = v1[n1];
533                 quad[2] = v1[p1];
534                 quad[3] = v2[p2];
535
536                 return checkValidQuad(vertex_storage, quad);
537         }
538         else if (v1[n1] == v2[p2]) {
539                 quad[0] = v1[current];
540                 quad[1] = v2[n2];
541                 quad[2] = v1[n1];
542                 quad[3] = v1[p1];
543
544                 return checkValidQuad(vertex_storage, quad);
545         }
546
547         return 0;
548 }
549
550 static bool Carve_checkDegeneratedFace(MeshSet<3>::face_t *face)
551 {
552         /* only tris and quads for now */
553         if (face->n_edges == 3) {
554                 return triangleArea(face->edge->prev->vert->v, face->edge->vert->v, face->edge->next->vert->v) < DBL_EPSILON;
555         }
556         else if (face->n_edges == 4) {
557                 return triangleArea(face->edge->vert->v, face->edge->next->vert->v, face->edge->next->next->vert->v) +
558                        triangleArea(face->edge->prev->vert->v, face->edge->vert->v, face->edge->next->next->vert->v) < DBL_EPSILON;
559         }
560
561         return false;
562 }
563
564 static BSP_CSGMesh *Carve_exportMesh(MeshSet<3>* &poly, carve::interpolate::FaceAttr<uint> &oface_num,
565                                      uint num_origfaces)
566 {
567         uint i;
568         BSP_CSGMesh *outputMesh = BSP_CSGMesh::New();
569
570         if (outputMesh == NULL)
571                 return NULL;
572
573         std::vector<BSP_MVertex> *vertices = new std::vector<BSP_MVertex>;
574
575         outputMesh->SetVertices(vertices);
576
577         std::map<MeshSet<3>::vertex_t*, uint> vertexToIndex_map;
578         std::vector<MeshSet<3>::vertex_t>::iterator it = poly->vertex_storage.begin();
579         for (i = 0; it != poly->vertex_storage.end(); ++i, ++it) {
580                 MeshSet<3>::vertex_t *vertex = &(*it);
581                 vertexToIndex_map[vertex] = i;
582         }
583
584         for (i = 0; i < poly->vertex_storage.size(); ++i ) {
585                 BSP_MVertex outVtx(MT_Point3 (poly->vertex_storage[i].v[0],
586                                               poly->vertex_storage[i].v[1],
587                                               poly->vertex_storage[i].v[2]));
588                 outVtx.m_edges.clear();
589                 outputMesh->VertexSet().push_back(outVtx);
590         }
591
592         // build vectors of faces for each original face and each vertex 
593         std::vector<std::vector<uint> > vi(poly->vertex_storage.size());
594         std::vector<std::vector<uint> > ofaces(num_origfaces);
595         MeshSet<3>::face_iter face_iter = poly->faceBegin();
596         for (i = 0; face_iter != poly->faceEnd(); ++face_iter, ++i) {
597                 MeshSet<3>::face_t *f = *face_iter;
598                 ofaces[oface_num.getAttribute(f)].push_back(i);
599                 MeshSet<3>::face_t::edge_iter_t edge_iter = f->begin();
600                 for (; edge_iter != f->end(); ++edge_iter) {
601                         int index = vertexToIndex_map[edge_iter->vert];
602                         vi[index].push_back(i);
603                 }
604         }
605
606         uint quadverts[4] = {0, 0, 0, 0};
607         // go over each set of faces which belong to an original face
608         std::vector< std::vector<uint> >::const_iterator fii;
609         uint orig = 0;
610         for (fii=ofaces.begin(); fii!=ofaces.end(); ++fii, ++orig) {
611                 std::vector<uint> fl = *fii;
612                 // go over a single set from an original face
613                 while (fl.size() > 0) {
614                         // remove one new face
615                         uint findex = fl.back();
616                         fl.pop_back();
617
618                         MeshSet<3>::face_t *f = *(poly->faceBegin() + findex);
619
620                         // for each vertex of this face, check other faces containing
621                         // that vertex to see if there is a neighbor also belonging to
622                         // the original face
623                         uint result = 0;
624
625                         MeshSet<3>::face_t::edge_iter_t edge_iter = f->begin();
626                         for (; edge_iter != f->end(); ++edge_iter) {
627                                 int v = vertexToIndex_map[edge_iter->vert];
628                                 for (uint pos2=0; !result && pos2 < vi[v].size();pos2++) {
629
630                                         // if we find the current face, ignore it
631                                         uint otherf = vi[v][pos2];
632                                         if (findex == otherf)
633                                                 continue;
634
635                                         MeshSet<3>::face_t *f2 = *(poly->faceBegin() + otherf);
636
637                                         // if other face doesn't have the same original face,
638                                         // ignore it also
639                                         uint other_orig = oface_num.getAttribute(f2);
640                                         if (orig != other_orig)
641                                                 continue;
642
643                                         // if, for some reason, we don't find the other face in
644                                         // the current set of faces, ignore it
645                                         uint other_index = 0;
646                                         while (other_index < fl.size() && fl[other_index] != otherf) ++other_index;
647                                         if (other_index == fl.size()) continue;
648
649                                         // see if the faces share an edge
650                                         result = quadMerge(&vertexToIndex_map, poly->vertex_storage, f, f2, v, quadverts);
651                                         // if faces can be merged, then remove the other face
652                                         // from the current set
653                                         if (result) {
654                                                 uint replace = fl.back();
655                                                 fl.pop_back();
656                                                 if(otherf != replace)
657                                                         fl[other_index] = replace;
658                                         }
659                                 }
660                         }
661
662                         bool degenerativeFace = false;
663
664                         if (!result) {
665                                 /* merged triangles are already checked for degenerative quad */
666                                 degenerativeFace = Carve_checkDegeneratedFace(f);
667                         }
668
669                         if (!degenerativeFace) {
670                                 // add all information except vertices to the output mesh
671                                 outputMesh->FaceSet().push_back(BSP_MFace());
672                                 BSP_MFace& outFace = outputMesh->FaceSet().back();
673                                 outFace.m_verts.clear();
674                                 outFace.m_plane.setValue(f->plane.N.v);
675                                 outFace.m_orig_face = orig;
676
677                                 // if we merged faces, use the list of common vertices; otherwise
678                                 // use the faces's vertices
679                                 if (result) {
680                                         // make quat using verts stored in result
681                                         outFace.m_verts.push_back(quadverts[0]);
682                                         outFace.m_verts.push_back(quadverts[1]);
683                                         outFace.m_verts.push_back(quadverts[2]);
684                                         outFace.m_verts.push_back(quadverts[3]);
685                                 } else {
686                                         MeshSet<3>::face_t::edge_iter_t edge_iter = f->begin();
687                                         for (; edge_iter != f->end(); ++edge_iter) {
688                                                 //int index = ofacevert_num.getAttribute(f, edge_iter.idx());
689                                                 int index = vertexToIndex_map[edge_iter->vert];
690                                                 outFace.m_verts.push_back( index );
691                                         }
692                                 }
693                         }
694                 }
695         }
696
697         // Build the mesh edges using topological informtion
698         outputMesh->BuildEdges();
699         
700         return outputMesh;
701 }
702
703 /**
704  * Performs a generic booleam operation, the entry point for external modules.
705  * @param opType Boolean operation type BOP_INTERSECTION, BOP_UNION, BOP_DIFFERENCE
706  * @param outputMesh Output mesh, the final result (the object C)
707  * @param obAFaces Object A faces list
708  * @param obAVertices Object A vertices list
709  * @param obBFaces Object B faces list
710  * @param obBVertices Object B vertices list
711  * @param interpFunc Interpolating function
712  * @return operation state: BOP_OK, BOP_NO_SOLID, BOP_ERROR
713  */
714 BoolOpState BOP_performBooleanOperation(BoolOpType                    opType,
715                                         BSP_CSGMesh**                 outputMesh,
716                                         CSG_FaceIteratorDescriptor    obAFaces,
717                                         CSG_VertexIteratorDescriptor  obAVertices,
718                                         CSG_FaceIteratorDescriptor    obBFaces,
719                                         CSG_VertexIteratorDescriptor  obBVertices)
720 {
721         carve::csg::CSG::OP op;
722         MeshSet<3> *left, *right, *output = NULL;
723         carve::csg::CSG csg;
724         carve::geom3d::Vector min, max;
725         carve::interpolate::FaceAttr<uint> oface_num;
726         uint num_origfaces = 0;
727
728         switch (opType) {
729                 case BOP_UNION:
730                         op = carve::csg::CSG::UNION;
731                         break;
732                 case BOP_INTERSECTION:
733                         op = carve::csg::CSG::INTERSECTION;
734                         break;
735                 case BOP_DIFFERENCE:
736                         op = carve::csg::CSG::A_MINUS_B;
737                         break;
738                 default:
739                         return BOP_ERROR;
740         }
741
742         left = Carve_addMesh(obAFaces, obAVertices, oface_num, num_origfaces );
743         right = Carve_addMesh(obBFaces, obBVertices, oface_num, num_origfaces );
744
745         min.x = max.x = left->vertex_storage[0].v.x;
746         min.y = max.y = left->vertex_storage[0].v.y;
747         min.z = max.z = left->vertex_storage[0].v.z;
748         for (uint i = 1; i < left->vertex_storage.size(); ++i) {
749                 min.x = MIN(min.x,left->vertex_storage[i].v.x);
750                 min.y = MIN(min.y,left->vertex_storage[i].v.y);
751                 min.z = MIN(min.z,left->vertex_storage[i].v.z);
752                 max.x = MAX(max.x,left->vertex_storage[i].v.x);
753                 max.y = MAX(max.y,left->vertex_storage[i].v.y);
754                 max.z = MAX(max.z,left->vertex_storage[i].v.z);
755         }
756         for (uint i = 0; i < right->vertex_storage.size(); ++i) {
757                 min.x = MIN(min.x,right->vertex_storage[i].v.x);
758                 min.y = MIN(min.y,right->vertex_storage[i].v.y);
759                 min.z = MIN(min.z,right->vertex_storage[i].v.z);
760                 max.x = MAX(max.x,right->vertex_storage[i].v.x);
761                 max.y = MAX(max.y,right->vertex_storage[i].v.y);
762                 max.z = MAX(max.z,right->vertex_storage[i].v.z);
763         }
764
765         carve::rescale::rescale scaler(min.x, min.y, min.z, max.x, max.y, max.z);
766         carve::rescale::fwd fwd_r(scaler);
767         carve::rescale::rev rev_r(scaler); 
768
769         left->transform(fwd_r);
770         right->transform(fwd_r);
771
772         // prepare operands for actual boolean operation. it's needed because operands might consist of
773         // several intersecting meshes and in case if another operands intersect an edge loop of intersecting that
774         // meshes tessellation of operation result can't be done properly. the only way to make such situations
775         // working is to union intersecting meshes of the same operand
776         Carve_unionIntersections(&left, &right, oface_num);
777
778         if(left->meshes.size() == 0 || right->meshes.size()==0) {
779                 // normally sohuldn't happen (zero-faces objects are handled by modifier itself), but
780                 // unioning intersecting meshes which doesn't have consistent normals might lead to
781                 // empty result which wouldn't work here
782
783                 delete left;
784                 delete right;
785
786                 return BOP_ERROR;
787         }
788
789         csg.hooks.registerHook(new carve::csg::CarveTriangulator, carve::csg::CSG::Hooks::PROCESS_OUTPUT_FACE_BIT);
790
791         oface_num.installHooks(csg);
792
793         try {
794                 output = csg.compute(left, right, op, NULL, carve::csg::CSG::CLASSIFY_EDGE);
795         }
796         catch(carve::exception e) {
797                 std::cerr << "CSG failed, exception " << e.str() << std::endl;
798         }
799         catch(...) {
800                 delete left;
801                 delete right;
802
803                 throw "Unknown error in Carve library";
804         }
805
806         delete left;
807         delete right;
808
809         if(!output)
810                 return BOP_ERROR;
811
812         output->transform(rev_r);
813
814         *outputMesh = Carve_exportMesh(output, oface_num, num_origfaces);
815         delete output;
816
817         return BOP_OK;
818 }