bb3a783548c93e3e962bd31fcc979c940ba77d3f
[blender.git] / intern / bsp / 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 bsp/intern/BOP_CarveInterface.cpp
30  *  \ingroup bsp
31  */
32
33 #include "BOP_Interface.h"
34 #include "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-5 * cross.length();
62
63         return fabsf(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(std::map<MeshSet<3>::vertex_t*, uint> *vertexToIndex_map, MeshSet<3>::face_t *face)
551 {
552         /* only tris and quads for now */
553         if (face->n_edges == 3) {
554                 uint v1, v2, v3;
555
556                 v1 = vertexToIndex_map->find(face->edge->prev->vert)->second;
557                 v2 = vertexToIndex_map->find(face->edge->vert)->second;
558                 v3 = vertexToIndex_map->find(face->edge->next->vert)->second;
559
560                 if (v1 == v2 || v2 == v3 || v1 == v3)
561                         return true;
562         }
563         else if (face->n_edges == 4) {
564                 uint v1, v2, v3, v4;
565
566                 v1 = vertexToIndex_map->find(face->edge->prev->vert)->second;
567                 v2 = vertexToIndex_map->find(face->edge->vert)->second;
568                 v3 = vertexToIndex_map->find(face->edge->next->vert)->second;
569                 v4 = vertexToIndex_map->find(face->edge->next->next->vert)->second;
570
571                 if (v1 == v2 || v1 == v3 || v1 == v4 || v2 == v3 || v2 == v4 || v3 == v4)
572                         return true;
573         }
574
575         return false;
576 }
577
578 static BSP_CSGMesh *Carve_exportMesh(MeshSet<3>* &poly, carve::interpolate::FaceAttr<uint> &oface_num,
579                                      uint num_origfaces)
580 {
581         uint i;
582         BSP_CSGMesh *outputMesh = BSP_CSGMesh::New();
583
584         if (outputMesh == NULL)
585                 return NULL;
586
587         std::vector<BSP_MVertex> *vertices = new std::vector<BSP_MVertex>;
588
589         outputMesh->SetVertices(vertices);
590
591         std::map<MeshSet<3>::vertex_t*, uint> vertexToIndex_map;
592         std::vector<MeshSet<3>::vertex_t>::iterator it = poly->vertex_storage.begin();
593         for (i = 0; it != poly->vertex_storage.end(); ++i, ++it) {
594                 MeshSet<3>::vertex_t *vertex = &(*it);
595                 vertexToIndex_map[vertex] = i;
596         }
597
598         for (i = 0; i < poly->vertex_storage.size(); ++i ) {
599                 BSP_MVertex outVtx(MT_Point3 (poly->vertex_storage[i].v[0],
600                                               poly->vertex_storage[i].v[1],
601                                               poly->vertex_storage[i].v[2]));
602                 outVtx.m_edges.clear();
603                 outputMesh->VertexSet().push_back(outVtx);
604         }
605
606         // build vectors of faces for each original face and each vertex 
607         std::vector<std::vector<uint> > vi(poly->vertex_storage.size());
608         std::vector<std::vector<uint> > ofaces(num_origfaces);
609         MeshSet<3>::face_iter face_iter = poly->faceBegin();
610         for (i = 0; face_iter != poly->faceEnd(); ++face_iter, ++i) {
611                 MeshSet<3>::face_t *f = *face_iter;
612
613                 if (Carve_checkDegeneratedFace(&vertexToIndex_map, f))
614                         continue;
615
616                 ofaces[oface_num.getAttribute(f)].push_back(i);
617
618                 MeshSet<3>::face_t::edge_iter_t edge_iter = f->begin();
619
620                 for (; edge_iter != f->end(); ++edge_iter) {
621                         int index = vertexToIndex_map[edge_iter->vert];
622                         vi[index].push_back(i);
623                 }
624         }
625
626         uint quadverts[4] = {0, 0, 0, 0};
627         // go over each set of faces which belong to an original face
628         std::vector< std::vector<uint> >::const_iterator fii;
629         uint orig = 0;
630         for (fii=ofaces.begin(); fii!=ofaces.end(); ++fii, ++orig) {
631                 std::vector<uint> fl = *fii;
632                 // go over a single set from an original face
633                 while (fl.size() > 0) {
634                         // remove one new face
635                         uint findex = fl.back();
636                         fl.pop_back();
637
638                         MeshSet<3>::face_t *f = *(poly->faceBegin() + findex);
639
640                         // for each vertex of this face, check other faces containing
641                         // that vertex to see if there is a neighbor also belonging to
642                         // the original face
643                         uint result = 0;
644
645                         MeshSet<3>::face_t::edge_iter_t edge_iter = f->begin();
646                         for (; edge_iter != f->end(); ++edge_iter) {
647                                 int v = vertexToIndex_map[edge_iter->vert];
648                                 for (uint pos2=0; !result && pos2 < vi[v].size();pos2++) {
649
650                                         // if we find the current face, ignore it
651                                         uint otherf = vi[v][pos2];
652                                         if (findex == otherf)
653                                                 continue;
654
655                                         MeshSet<3>::face_t *f2 = *(poly->faceBegin() + otherf);
656
657                                         // if other face doesn't have the same original face,
658                                         // ignore it also
659                                         uint other_orig = oface_num.getAttribute(f2);
660                                         if (orig != other_orig)
661                                                 continue;
662
663                                         // if, for some reason, we don't find the other face in
664                                         // the current set of faces, ignore it
665                                         uint other_index = 0;
666                                         while (other_index < fl.size() && fl[other_index] != otherf) ++other_index;
667                                         if (other_index == fl.size()) continue;
668
669                                         // see if the faces share an edge
670                                         result = quadMerge(&vertexToIndex_map, poly->vertex_storage, f, f2, v, quadverts);
671                                         // if faces can be merged, then remove the other face
672                                         // from the current set
673                                         if (result) {
674                                                 uint replace = fl.back();
675                                                 fl.pop_back();
676                                                 if(otherf != replace)
677                                                         fl[other_index] = replace;
678                                         }
679                                 }
680                         }
681
682                         // add all information except vertices to the output mesh
683                         outputMesh->FaceSet().push_back(BSP_MFace());
684                         BSP_MFace& outFace = outputMesh->FaceSet().back();
685                         outFace.m_verts.clear();
686                         outFace.m_plane.setValue(f->plane.N.v);
687                         outFace.m_orig_face = orig;
688
689                         // if we merged faces, use the list of common vertices; otherwise
690                         // use the faces's vertices
691                         if (result) {
692                                 // make quat using verts stored in result
693                                 outFace.m_verts.push_back(quadverts[0]);
694                                 outFace.m_verts.push_back(quadverts[1]);
695                                 outFace.m_verts.push_back(quadverts[2]);
696                                 outFace.m_verts.push_back(quadverts[3]);
697                         } else {
698                                 MeshSet<3>::face_t::edge_iter_t edge_iter = f->begin();
699                                 for (; edge_iter != f->end(); ++edge_iter) {
700                                         //int index = ofacevert_num.getAttribute(f, edge_iter.idx());
701                                         int index = vertexToIndex_map[edge_iter->vert];
702                                         outFace.m_verts.push_back( index );
703                                 }
704                         }
705                 }
706         }
707
708         // Build the mesh edges using topological informtion
709         outputMesh->BuildEdges();
710         
711         return outputMesh;
712 }
713
714 /**
715  * Performs a generic booleam operation, the entry point for external modules.
716  * @param opType Boolean operation type BOP_INTERSECTION, BOP_UNION, BOP_DIFFERENCE
717  * @param outputMesh Output mesh, the final result (the object C)
718  * @param obAFaces Object A faces list
719  * @param obAVertices Object A vertices list
720  * @param obBFaces Object B faces list
721  * @param obBVertices Object B vertices list
722  * @param interpFunc Interpolating function
723  * @return operation state: BOP_OK, BOP_NO_SOLID, BOP_ERROR
724  */
725 BoolOpState BOP_performBooleanOperation(BoolOpType                    opType,
726                                         BSP_CSGMesh**                 outputMesh,
727                                         CSG_FaceIteratorDescriptor    obAFaces,
728                                         CSG_VertexIteratorDescriptor  obAVertices,
729                                         CSG_FaceIteratorDescriptor    obBFaces,
730                                         CSG_VertexIteratorDescriptor  obBVertices)
731 {
732         carve::csg::CSG::OP op;
733         MeshSet<3> *left, *right, *output = NULL;
734         carve::csg::CSG csg;
735         carve::geom3d::Vector min, max;
736         carve::interpolate::FaceAttr<uint> oface_num;
737         uint num_origfaces = 0;
738
739         switch (opType) {
740                 case BOP_UNION:
741                         op = carve::csg::CSG::UNION;
742                         break;
743                 case BOP_INTERSECTION:
744                         op = carve::csg::CSG::INTERSECTION;
745                         break;
746                 case BOP_DIFFERENCE:
747                         op = carve::csg::CSG::A_MINUS_B;
748                         break;
749                 default:
750                         return BOP_ERROR;
751         }
752
753         left = Carve_addMesh(obAFaces, obAVertices, oface_num, num_origfaces );
754         right = Carve_addMesh(obBFaces, obBVertices, oface_num, num_origfaces );
755
756         min.x = max.x = left->vertex_storage[0].v.x;
757         min.y = max.y = left->vertex_storage[0].v.y;
758         min.z = max.z = left->vertex_storage[0].v.z;
759         for (uint i = 1; i < left->vertex_storage.size(); ++i) {
760                 min.x = MIN(min.x,left->vertex_storage[i].v.x);
761                 min.y = MIN(min.y,left->vertex_storage[i].v.y);
762                 min.z = MIN(min.z,left->vertex_storage[i].v.z);
763                 max.x = MAX(max.x,left->vertex_storage[i].v.x);
764                 max.y = MAX(max.y,left->vertex_storage[i].v.y);
765                 max.z = MAX(max.z,left->vertex_storage[i].v.z);
766         }
767         for (uint i = 0; i < right->vertex_storage.size(); ++i) {
768                 min.x = MIN(min.x,right->vertex_storage[i].v.x);
769                 min.y = MIN(min.y,right->vertex_storage[i].v.y);
770                 min.z = MIN(min.z,right->vertex_storage[i].v.z);
771                 max.x = MAX(max.x,right->vertex_storage[i].v.x);
772                 max.y = MAX(max.y,right->vertex_storage[i].v.y);
773                 max.z = MAX(max.z,right->vertex_storage[i].v.z);
774         }
775
776         carve::rescale::rescale scaler(min.x, min.y, min.z, max.x, max.y, max.z);
777         carve::rescale::fwd fwd_r(scaler);
778         carve::rescale::rev rev_r(scaler); 
779
780         left->transform(fwd_r);
781         right->transform(fwd_r);
782
783         // prepare operands for actual boolean operation. it's needed because operands might consist of
784         // several intersecting meshes and in case if another operands intersect an edge loop of intersecting that
785         // meshes tessellation of operation result can't be done properly. the only way to make such situations
786         // working is to union intersecting meshes of the same operand
787         Carve_unionIntersections(&left, &right, oface_num);
788
789         if(left->meshes.size() == 0 || right->meshes.size()==0) {
790                 // normally sohuldn't happen (zero-faces objects are handled by modifier itself), but
791                 // unioning intersecting meshes which doesn't have consistent normals might lead to
792                 // empty result which wouldn't work here
793
794                 delete left;
795                 delete right;
796
797                 return BOP_ERROR;
798         }
799
800         csg.hooks.registerHook(new carve::csg::CarveTriangulator, carve::csg::CSG::Hooks::PROCESS_OUTPUT_FACE_BIT);
801
802         oface_num.installHooks(csg);
803
804         try {
805                 output = csg.compute(left, right, op, NULL, carve::csg::CSG::CLASSIFY_EDGE);
806         }
807         catch(carve::exception e) {
808                 std::cerr << "CSG failed, exception " << e.str() << std::endl;
809         }
810         catch(...) {
811                 delete left;
812                 delete right;
813
814                 throw "Unknown error in Carve library";
815         }
816
817         delete left;
818         delete right;
819
820         if(!output)
821                 return BOP_ERROR;
822
823         output->transform(rev_r);
824
825         *outputMesh = Carve_exportMesh(output, oface_num, num_origfaces);
826         delete output;
827
828         return BOP_OK;
829 }