Rework carve integration into boolean modifier

Goal of this commit is to support NGons for boolean modifier
(currently mesh is being tessellated before performing boolean
operation) and also solve the limitation of loosing edge custom
data layers after boolean operation is performed.

Main idea is to make it so boolean modifier uses Carve library
directly via it's C-API, avoiding BSP intermediate level which
was doubling amount of memory needed for the operation and which
also used quite reasonable amount of overhead time.

Perhaps memory usage and CPU usage are the same after all the
features are implemented but we've got support now:

- ORIGINDEX for all the geometry
- Interpolation of edge custom data (seams, crease)
- NGons support

Triangulation rule is changed now as well, so now non-flat
polygons are not being merged back after Carve work. This is
so because it's not so trivial to support for NGons and
having different behavior for quads and NGons is even more
creepy.

Reviewers: lukastoenne, campbellbarton

Differential Revision: https://developer.blender.org/D274

9 files changed, 1666 insertions, 1 deletions

diff --git a/extern/carve/CMakeLists.txt b/extern/carve/CMakeLists.txt
index 5e917ac1e44..11a92f92abc 100644
--- a/extern/carve/CMakeLists.txt
+++ b/extern/carve/CMakeLists.txt
@@ -35,6 +35,8 @@ set(INC_SYS
 )
 
 set(SRC
+	carve-capi.cc
+	carve-util.cc
 	lib/aabb.cpp
 	lib/carve.cpp
 	lib/convex_hull.cpp
@@ -62,6 +64,8 @@ set(SRC
 	lib/timing.cpp
 	lib/triangulator.cpp
 
+	carve-capi.h
+	carve-util.h
 	lib/csg_collector.hpp
 	lib/csg_data.hpp
 	lib/csg_detail.hpp
diff --git a/extern/carve/bundle.sh b/extern/carve/bundle.sh
index 05967d6131d..91d5f44a880 100755
--- a/extern/carve/bundle.sh
+++ b/extern/carve/bundle.sh
@@ -70,8 +70,12 @@ set(INC_SYS
 )
 
 set(SRC
+	carve-capi.cc
+	carve-util.cc
 ${sources}
 
+	carve-capi.h
+	carve-util.h
 ${headers}
 
 ${includes}
diff --git a/extern/carve/carve-capi.cc b/extern/carve/carve-capi.cc
new file mode 100644
index 00000000000..7478c341663
--- /dev/null
+++ b/extern/carve/carve-capi.cc
@@ -0,0 +1,580 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2014 Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Blender Foundation,
+ *                 Sergey Sharybin
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#include "carve-capi.h"
+#include "carve-util.h"
+
+#include <carve/interpolator.hpp>
+#include <carve/rescale.hpp>
+
+using carve::mesh::MeshSet;
+
+typedef std::pair<int, int> OrigIndex;
+typedef std::pair<MeshSet<3>::vertex_t *, MeshSet<3>::vertex_t *> VertexPair;
+typedef carve::interpolate::VertexAttr<OrigIndex> OrigVertMapping;
+typedef carve::interpolate::FaceAttr<OrigIndex> OrigFaceMapping;
+typedef carve::interpolate::FaceEdgeAttr<OrigIndex> OrigFaceEdgeMapping;
+
+typedef struct CarveMeshDescr {
+	// Stores mesh data itself.
+	MeshSet<3> *poly;
+
+	// N-th element of the vector indicates index of an original mesh loop.
+	std::vector<int> orig_loop_index_map;
+
+	// N-th element of the vector indicates index of an original mesh poly.
+	std::vector<int> orig_poly_index_map;
+
+	// The folloving mapping is only filled in for output mesh.
+
+	// Mapping from the face verts back to original vert index.
+	OrigVertMapping orig_vert_mapping;
+
+	// Mapping from the face edges back to (original edge index, original loop index).
+	OrigFaceEdgeMapping orig_face_edge_mapping;
+
+	// Mapping from the faces back to original poly index.
+	OrigFaceMapping orig_face_mapping;
+} CarveMeshDescr;
+
+namespace {
+
+template <typename T1, typename T2>
+void edgeIndexMap_put(std::unordered_map<std::pair<T1, T1>, T2> *edge_map,
+                      const T1 &v1,
+                      const T1 &v2,
+                      const T2 &index)
+{
+	if (v1 < v2) {
+		(*edge_map)[std::make_pair(v1, v2)] = index;
+	}
+	else {
+		(*edge_map)[std::make_pair(v2, v1)] = index;
+	}
+}
+
+template <typename T1, typename T2>
+const T2 &edgeIndexMap_get(const std::unordered_map<std::pair<T1, T1>, T2> &edge_map,
+                           const T1 &v1,
+                           const T1 &v2)
+{
+	typedef std::unordered_map<std::pair<T1, T1>, T2> Map;
+	typename Map::const_iterator found;
+
+	if (v1 < v2) {
+		found = edge_map.find(std::make_pair(v1, v2));
+	}
+	else {
+		found = edge_map.find(std::make_pair(v2, v1));
+	}
+
+	assert(found != edge_map.end());
+	return found->second;
+}
+
+template <typename T1, typename T2>
+bool edgeIndexMap_get_if_exists(const std::unordered_map<std::pair<T1, T1>, T2> &edge_map,
+                                const T1 &v1,
+                                const T1 &v2,
+                                T2 *out)
+{
+	typedef std::unordered_map<std::pair<T1, T1>, T2> Map;
+	typename Map::const_iterator found;
+
+	if (v1 < v2) {
+		found = edge_map.find(std::make_pair(v1, v2));
+	}
+	else {
+		found = edge_map.find(std::make_pair(v2, v1));
+	}
+
+	if (found == edge_map.end()) {
+		return false;
+	}
+	*out = found->second;
+	return true;
+}
+
+template <typename T>
+inline int indexOf(const T *element, const std::vector<T> &vector_from)
+{
+	return element - &vector_from.at(0);
+}
+
+void initOrigIndexMeshFaceMapping(CarveMeshDescr *mesh,
+                                  int which_mesh,
+                                  const std::vector<int> &orig_loop_index_map,
+                                  const std::vector<int> &orig_poly_index_map,
+                                  OrigVertMapping *orig_vert_mapping,
+                                  OrigFaceEdgeMapping *orig_face_edge_mapping,
+                                  OrigFaceMapping *orig_face_attr)
+{
+	MeshSet<3> *poly = mesh->poly;
+
+	std::vector<MeshSet<3>::vertex_t>::iterator vertex_iter =
+		poly->vertex_storage.begin();
+	for (int i = 0;
+	     vertex_iter != poly->vertex_storage.end();
+	     ++i, ++vertex_iter)
+	{
+		MeshSet<3>::vertex_t *vertex = &(*vertex_iter);
+		orig_vert_mapping->setAttribute(vertex,
+		                                std::make_pair(which_mesh, i));
+	}
+
+	MeshSet<3>::face_iter face_iter = poly->faceBegin();
+	for (int i = 0, loop_map_index = 0;
+	     face_iter != poly->faceEnd();
+	     ++face_iter, ++i)
+	{
+		const MeshSet<3>::face_t *face = *face_iter;
+
+		// Mapping from carve face back to original poly index.
+		int orig_poly_index = orig_poly_index_map[i];
+		orig_face_attr->setAttribute(face, std::make_pair(which_mesh, orig_poly_index));
+
+		for (MeshSet<3>::face_t::const_edge_iter_t edge_iter = face->begin();
+		     edge_iter != face->end();
+		     ++edge_iter, ++loop_map_index)
+		{
+			int orig_loop_index = orig_loop_index_map[loop_map_index];
+			if (orig_loop_index != -1) {
+				// Mapping from carve face edge back to original loop index.
+				orig_face_edge_mapping->setAttribute(face,
+				                                     edge_iter.idx(),
+				                                     std::make_pair(which_mesh,
+				                                                    orig_loop_index));
+			}
+		}
+	}
+}
+
+void initOrigIndexMapping(CarveMeshDescr *left_mesh,
+                          CarveMeshDescr *right_mesh,
+                          OrigVertMapping *orig_vert_mapping,
+                          OrigFaceEdgeMapping *orig_face_edge_mapping,
+                          OrigFaceMapping *orig_face_mapping)
+{
+	initOrigIndexMeshFaceMapping(left_mesh,
+	                             CARVE_MESH_LEFT,
+	                             left_mesh->orig_loop_index_map,
+	                             left_mesh->orig_poly_index_map,
+	                             orig_vert_mapping,
+	                             orig_face_edge_mapping,
+	                             orig_face_mapping);
+
+	initOrigIndexMeshFaceMapping(right_mesh,
+	                             CARVE_MESH_RIGHT,
+	                             right_mesh->orig_loop_index_map,
+	                             right_mesh->orig_poly_index_map,
+	                             orig_vert_mapping,
+	                             orig_face_edge_mapping,
+	                             orig_face_mapping);
+}
+
+}  // namespace
+
+CarveMeshDescr *carve_addMesh(struct ImportMeshData *import_data,
+                              CarveMeshImporter *mesh_importer)
+{
+#define MAX_STATIC_VERTS 64
+
+	CarveMeshDescr *mesh_descr = new CarveMeshDescr;
+
+	// Import verices from external mesh to Carve.
+	int num_verts = mesh_importer->getNumVerts(import_data);
+	std::vector<carve::geom3d::Vector> vertices;
+	vertices.reserve(num_verts);
+	for (int i = 0; i < num_verts; i++) {
+		float position[3];
+		mesh_importer->getVertCoord(import_data, i, position);
+		vertices.push_back(carve::geom::VECTOR(position[0],
+		                                       position[1],
+		                                       position[2]));
+	}
+
+	// Import polys from external mesh to Carve.
+	int verts_of_poly_static[MAX_STATIC_VERTS];
+	int *verts_of_poly_dynamic = NULL;
+	int verts_of_poly_dynamic_size = 0;
+
+	int num_loops = mesh_importer->getNumLoops(import_data);
+	int num_polys = mesh_importer->getNumPolys(import_data);
+	int loop_index = 0;
+	int num_tessellated_polys = 0;
+	std::vector<int> face_indices;
+	face_indices.reserve(num_loops);
+	mesh_descr->orig_loop_index_map.reserve(num_polys);
+	mesh_descr->orig_poly_index_map.reserve(num_polys);
+	for (int i = 0; i < num_polys; i++) {
+		int verts_per_poly =
+			mesh_importer->getNumPolyVerts(import_data, i);
+		int *verts_of_poly;
+
+		if (verts_per_poly <= MAX_STATIC_VERTS) {
+			verts_of_poly = verts_of_poly_static;
+		}
+		else {
+			if (verts_of_poly_dynamic_size < verts_per_poly) {
+				if (verts_of_poly_dynamic != NULL) {
+					delete [] verts_of_poly_dynamic;
+				}
+				verts_of_poly_dynamic = new int[verts_per_poly];
+				verts_of_poly_dynamic_size = verts_per_poly;
+			}
+			verts_of_poly = verts_of_poly_dynamic;
+		}
+
+		mesh_importer->getPolyVerts(import_data, i, verts_of_poly);
+
+		carve::math::Matrix3 axis_matrix;
+		if (!carve_checkPolyPlanarAndGetNormal(vertices,
+		                                       verts_per_poly,
+		                                       verts_of_poly,
+		                                       &axis_matrix)) {
+			int num_triangles;
+
+			num_triangles = carve_triangulatePoly(import_data,
+			                                      mesh_importer,
+			                                      i,
+			                                      loop_index,
+			                                      vertices,
+			                                      verts_per_poly,
+			                                      verts_of_poly,
+			                                      axis_matrix,
+			                                      &face_indices,
+			                                      &mesh_descr->orig_loop_index_map,
+			                                      &mesh_descr->orig_poly_index_map);
+
+			num_tessellated_polys += num_triangles;
+		}
+		else {
+			face_indices.push_back(verts_per_poly);
+			for (int j = 0; j < verts_per_poly; ++j) {
+				mesh_descr->orig_loop_index_map.push_back(loop_index++);
+				face_indices.push_back(verts_of_poly[j]);
+			}
+			mesh_descr->orig_poly_index_map.push_back(i);
+			num_tessellated_polys++;
+		}
+	}
+
+	if (verts_of_poly_dynamic != NULL) {
+		delete [] verts_of_poly_dynamic;
+	}
+
+	mesh_descr->poly = new MeshSet<3> (vertices,
+	                                   num_tessellated_polys,
+	                                   face_indices);
+
+	return mesh_descr;
+
+#undef MAX_STATIC_VERTS
+}
+
+void carve_deleteMesh(CarveMeshDescr *mesh_descr)
+{
+	delete mesh_descr->poly;
+	delete mesh_descr;
+}
+
+bool carve_performBooleanOperation(CarveMeshDescr *left_mesh,
+                                   CarveMeshDescr *right_mesh,
+                                   int operation,
+                                   CarveMeshDescr **output_mesh)
+{
+	*output_mesh = NULL;
+
+	carve::csg::CSG::OP op;
+	switch (operation) {
+#define OP_CONVERT(the_op) \
+		case CARVE_OP_ ## the_op: \
+			op = carve::csg::CSG::the_op; \
+			break;
+		OP_CONVERT(UNION)
+		OP_CONVERT(INTERSECTION)
+		OP_CONVERT(A_MINUS_B)
+		default:
+			return false;
+#undef OP_CONVERT
+	}
+
+	CarveMeshDescr *output_descr = new CarveMeshDescr;
+	output_descr->poly = NULL;
+	try {
+		MeshSet<3> *left = left_mesh->poly, *right = right_mesh->poly;
+		carve::geom3d::Vector min, max;
+
+		// TODO(sergey): Make importer/exporter to care about re-scale
+		// to save extra mesh iteration here.
+		carve_getRescaleMinMax(left, right, &min, &max);
+
+		carve::rescale::rescale scaler(min.x, min.y, min.z, max.x, max.y, max.z);
+		carve::rescale::fwd fwd_r(scaler);
+		carve::rescale::rev rev_r(scaler);
+
+		left->transform(fwd_r);
+		right->transform(fwd_r);
+
+		// Initialize attributes for maping from boolean result mesh back to
+		// original geometry indices.
+		initOrigIndexMapping(left_mesh, right_mesh,
+		                     &output_descr->orig_vert_mapping,
+		                     &output_descr->orig_face_edge_mapping,
+		                     &output_descr->orig_face_mapping);
+
+		carve::csg::CSG csg;
+
+		output_descr->orig_vert_mapping.installHooks(csg);
+		output_descr->orig_face_edge_mapping.installHooks(csg);
+		output_descr->orig_face_mapping.installHooks(csg);
+
+		// Prepare operands for actual boolean operation.
+		//
+		// It's needed because operands might consist of several intersecting
+		// meshes and in case of another operands intersect an edge loop of
+		// intersecting that meshes tessellation of operation result can't be
+		// done properly. The only way to make such situations working is to
+		// union intersecting meshes of the same operand.
+		carve_unionIntersections(&csg, &left, &right);
+		left_mesh->poly = left;
+		right_mesh->poly = right;
+
+		if (left->meshes.size() == 0 || right->meshes.size() == 0) {
+			// Normally shouldn't happen (zero-faces objects are handled by
+			// modifier itself), but unioning intersecting meshes which doesn't
+			// have consistent normals might lead to empty result which
+			// wouldn't work here.
+
+			return false;
+		}
+
+		output_descr->poly = csg.compute(left,
+		                                 right,
+		                                 op,
+		                                 NULL,
+		                                 carve::csg::CSG::CLASSIFY_EDGE);
+		if (output_descr->poly) {
+			output_descr->poly->transform(rev_r);
+		}
+	}
+	catch (carve::exception e) {
+		std::cerr << "CSG failed, exception " << e.str() << std::endl;
+	}
+	catch (...) {
+		std::cerr << "Unknown error in Carve library" << std::endl;
+	}
+
+	*output_mesh = output_descr;
+
+	return output_descr->poly != NULL;
+}
+
+static void exportMesh_handle_edges_list(MeshSet<3> *poly,
+                                         const std::vector<MeshSet<3>::edge_t*> &edges,
+                                         int start_edge_index,
+                                         CarveMeshExporter *mesh_exporter,
+                                         struct ExportMeshData *export_data,
+                                         const std::unordered_map<VertexPair, OrigIndex> &edge_origindex_map,
+                                         std::unordered_map<VertexPair, int> *edge_map)
+{
+	for (int i = 0, edge_index = start_edge_index;
+	     i < edges.size();
+	     ++i, ++edge_index)
+	{
+		MeshSet<3>::edge_t *edge = edges.at(i);
+		MeshSet<3>::vertex_t *v1 = edge->vert;
+		MeshSet<3>::vertex_t *v2 = edge->next->vert;
+
+		OrigIndex orig_edge_index;
+
+		if (!edgeIndexMap_get_if_exists(edge_origindex_map,
+		                                v1,
+		                                v2,
+		                                &orig_edge_index))
+		{
+			orig_edge_index.first = CARVE_MESH_NONE;
+			orig_edge_index.second = -1;
+		}
+
+		mesh_exporter->setEdge(export_data,
+		                       edge_index,
+		                       indexOf(v1, poly->vertex_storage),
+		                       indexOf(v2, poly->vertex_storage),
+		                       orig_edge_index.first,
+		                       orig_edge_index.second);
+
+		edgeIndexMap_put(edge_map, v1, v2, edge_index);
+	}
+}
+
+void carve_exportMesh(CarveMeshDescr *mesh_descr,
+                      CarveMeshExporter *mesh_exporter,
+                      struct ExportMeshData *export_data)
+{
+	OrigIndex origindex_none = std::make_pair((int)CARVE_MESH_NONE, -1);
+	MeshSet<3> *poly = mesh_descr->poly;
+	int num_vertices = poly->vertex_storage.size();
+	int num_edges = 0, num_loops = 0, num_polys = 0;
+
+	// Count edges from all manifolds.
+	for (int i = 0; i < poly->meshes.size(); ++i) {
+		carve::mesh::Mesh<3> *mesh = poly->meshes[i];
+		num_edges += mesh->closed_edges.size() + mesh->open_edges.size();
+	}
+
+	// Count polys and loops from all manifolds.
+	for (MeshSet<3>::face_iter face_iter = poly->faceBegin();
+	     face_iter != poly->faceEnd();
+	     ++face_iter, ++num_polys)
+	{
+		MeshSet<3>::face_t *face = *face_iter;
+		num_loops += face->n_edges;
+	}
+
+	// Initialize arrays for geometry in exported mesh.
+	mesh_exporter->initGeomArrays(export_data,
+	                              num_vertices,
+	                              num_edges,
+	                              num_loops,
+	                              num_polys);
+
+	// Get mapping from edge denoted by vertex pair to original edge index,
+	//
+	// This is needed because internally Carve interpolates data for per-face
+	// edges rather then having some global edge storage.
+	std::unordered_map<VertexPair, OrigIndex> edge_origindex_map;
+	for (MeshSet<3>::face_iter face_iter = poly->faceBegin();
+	     face_iter != poly->faceEnd();
+	     ++face_iter)
+	{
+		MeshSet<3>::face_t *face = *face_iter;
+		for (MeshSet<3>::face_t::edge_iter_t edge_iter = face->begin();
+		     edge_iter != face->end();
+		     ++edge_iter)
+		{
+			MeshSet<3>::edge_t &edge = *edge_iter;
+			int edge_iter_index = edge_iter.idx();
+
+			const OrigIndex &orig_loop_index =
+				mesh_descr->orig_face_edge_mapping.getAttribute(face,
+				                                                edge_iter_index,
+				                                                origindex_none);
+
+			if (orig_loop_index.first != CARVE_MESH_NONE) {
+				OrigIndex orig_edge_index;
+
+				orig_edge_index.first = orig_loop_index.first;
+				orig_edge_index.second =
+					mesh_exporter->mapLoopToEdge(export_data,
+					                             orig_loop_index.first,
+					                             orig_loop_index.second);
+
+				MeshSet<3>::vertex_t *v1 = edge.vert;
+				MeshSet<3>::vertex_t *v2 = edge.next->vert;
+
+				edgeIndexMap_put(&edge_origindex_map, v1, v2, orig_edge_index);
+			}
+		}
+	}
+
+	// Export all the vertices.
+	std::vector<MeshSet<3>::vertex_t>::iterator vertex_iter = poly->vertex_storage.begin();
+	for (int i = 0; vertex_iter != poly->vertex_storage.end(); ++i, ++vertex_iter) {
+		MeshSet<3>::vertex_t *vertex = &(*vertex_iter);
+
+		OrigIndex orig_vert_index =
+			mesh_descr->orig_vert_mapping.getAttribute(vertex, origindex_none);
+
+		float coord[3];
+		coord[0] = vertex->v[0];
+		coord[1] = vertex->v[1];
+		coord[2] = vertex->v[2];
+		mesh_exporter->setVert(export_data, i, coord,
+		                       orig_vert_index.first,
+		                       orig_vert_index.second);
+	}
+
+	// Export all the edges.
+	std::unordered_map<VertexPair, int> edge_map;
+	for (int i = 0, edge_index = 0; i < poly->meshes.size(); ++i) {
+		carve::mesh::Mesh<3> *mesh = poly->meshes[i];
+
+		// Export closed edges.
+		exportMesh_handle_edges_list(poly,
+		                             mesh->closed_edges,
+		                             edge_index,
+		                             mesh_exporter,
+		                             export_data,
+		                             edge_origindex_map,
+		                             &edge_map);
+
+		// Export open edges.
+		exportMesh_handle_edges_list(poly,
+		                             mesh->open_edges,
+		                             edge_index,
+		                             mesh_exporter,
+		                             export_data,
+		                             edge_origindex_map,
+		                             &edge_map);
+
+		edge_index += mesh->closed_edges.size();
+	}
+
+	// Export all the loops and polys.
+	MeshSet<3>::face_iter face_iter = poly->faceBegin();
+	for (int loop_index = 0, poly_index = 0;
+	     face_iter != poly->faceEnd();
+	     ++face_iter, ++poly_index)
+	{
+		int start_loop_index = loop_index;
+		MeshSet<3>::face_t *face = *face_iter;
+		const OrigIndex &orig_face_index =
+			mesh_descr->orig_face_mapping.getAttribute(face, origindex_none);
+
+		for (MeshSet<3>::face_t::edge_iter_t edge_iter = face->begin();
+		     edge_iter != face->end();
+		     ++edge_iter, ++loop_index)
+		{
+			MeshSet<3>::edge_t &edge = *edge_iter;
+			const OrigIndex &orig_loop_index =
+				mesh_descr->orig_face_edge_mapping.getAttribute(face,
+				                                                edge_iter.idx(),
+				                                                origindex_none);
+
+			mesh_exporter->setLoop(export_data,
+		                           loop_index,
+		                           indexOf(edge.vert, poly->vertex_storage),
+		                           edgeIndexMap_get(edge_map, edge.vert, edge.next->vert),
+		                           orig_loop_index.first,
+		                           orig_loop_index.second);
+		}
+
+		mesh_exporter->setPoly(export_data,
+		                       poly_index, start_loop_index, face->n_edges,
+		                       orig_face_index.first, orig_face_index.second);
+	}
+}
diff --git a/extern/carve/carve-capi.h b/extern/carve/carve-capi.h
new file mode 100644
index 00000000000..25704dfeb48
--- /dev/null
+++ b/extern/carve/carve-capi.h
@@ -0,0 +1,164 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2014 Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Blender Foundation,
+ *                 Sergey Sharybin
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#ifndef __CARVE_CAPI_H__
+#define __CARVE_CAPI_H__
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct CarveMeshDescr;
+
+//
+// Importer from external storage to Carve module
+//
+
+struct ImportMeshData;
+
+// Get number of vertices.
+typedef int (*CarveImporter_GetNumVerts) (struct ImportMeshData *import_data);
+
+// Get number of edges.
+typedef int (*CarveImporter_GetNumEdges) (struct ImportMeshData *import_data);
+
+// Get number of loops.
+typedef int (*CarveImporter_GetNumLoops) (struct ImportMeshData *import_data);
+
+// Get number of polys.
+typedef int (*CarveImporter_GetNumPolys) (struct ImportMeshData *import_data);
+
+// Get 3D coordinate of vertex with given index.
+typedef void (*CarveImporter_GetVertCoord) (struct ImportMeshData *import_data, int vert_index, float coord[3]);
+
+// Get index of vertices which are adjucent to edge specified by it's index.
+typedef void (*CarveImporter_GetEdgeVerts) (struct ImportMeshData *import_data, int edge_index, int *v1, int *v2);
+
+// Get number of adjucent vertices to the poly specified by it's index.
+typedef int (*CarveImporter_GetPolyNumVerts) (struct ImportMeshData *import_data, int poly_index);
+
+// Get list of adjucent vertices to the poly specified by it's index.
+typedef void (*CarveImporter_GetPolyVerts) (struct ImportMeshData *import_data, int poly_index, int *verts);
+
+// Triangulate 2D polygon.
+typedef int (*CarveImporter_Triangulate2DPoly) (struct ImportMeshData *import_data,
+                                                const float (*vertices)[2], int num_vertices,
+                                                unsigned int (*triangles)[3]);
+
+typedef struct CarveMeshImporter {
+	CarveImporter_GetNumVerts getNumVerts;
+	CarveImporter_GetNumEdges getNumEdges;
+	CarveImporter_GetNumLoops getNumLoops;
+	CarveImporter_GetNumPolys getNumPolys;
+	CarveImporter_GetVertCoord getVertCoord;
+	CarveImporter_GetEdgeVerts getEdgeVerts;
+	CarveImporter_GetPolyNumVerts getNumPolyVerts;
+	CarveImporter_GetPolyVerts getPolyVerts;
+	CarveImporter_Triangulate2DPoly triangulate2DPoly;
+} CarveMeshImporter;
+
+//
+// Exporter from Carve module to external storage
+//
+
+struct ExportMeshData;
+
+// Initialize arrays for geometry.
+typedef void (*CarveExporter_InitGeomArrays) (struct ExportMeshData *export_data,
+                                              int num_verts, int num_edges,
+                                              int num_polys, int num_loops);
+
+// Set coordinate of vertex with given index.
+typedef void (*CarveExporter_SetVert) (struct ExportMeshData *export_data,
+                                       int vert_index, float coord[3],
+                                       int which_orig_mesh, int orig_edge_index);
+
+// Set vertices which are adjucent to the edge specified by it's index.
+typedef void (*CarveExporter_SetEdge) (struct ExportMeshData *export_data,
+                                       int edge_index, int v1, int v2,
+                                       int which_orig_mesh, int orig_edge_index);
+
+// Set adjucent loops to the poly specified by it's index.
+typedef void (*CarveExporter_SetPoly) (struct ExportMeshData *export_data,
+                                       int poly_index, int start_loop, int num_loops,
+                                       int which_orig_mesh, int orig_poly_index);
+
+// Set list vertex and edge which are adjucent to loop with given index.
+typedef void (*CarveExporter_SetLoop) (struct ExportMeshData *export_data,
+                                       int loop_index, int vertex, int edge,
+                                       int which_orig_mesh, int orig_loop_index);
+
+// Get edge index from a loop index for a given original mesh.
+//
+// A bit of a bummer to access original operands data on export stage,
+// but Blender side still does have this information in derived meshes
+// and we use API to get this data instead of duplicating it in Carve
+// API side. This is because of optimizations reasons.
+typedef int (*CarveExporter_MapLoopToEdge) (struct ExportMeshData *export_data,
+                                            int which_mesh, int loop_index);
+
+typedef struct CarveMeshExporter {
+	CarveExporter_InitGeomArrays initGeomArrays;
+	CarveExporter_SetVert setVert;
+	CarveExporter_SetEdge setEdge;
+	CarveExporter_SetPoly setPoly;
+	CarveExporter_SetLoop setLoop;
+	CarveExporter_MapLoopToEdge mapLoopToEdge;
+} CarveMeshExporter;
+
+enum {
+	CARVE_OP_UNION,
+	CARVE_OP_INTERSECTION,
+	CARVE_OP_A_MINUS_B,
+};
+
+enum {
+	CARVE_MESH_NONE,
+	CARVE_MESH_LEFT,
+	CARVE_MESH_RIGHT
+};
+
+struct CarveMeshDescr *carve_addMesh(struct ImportMeshData *import_data,
+                                     CarveMeshImporter *mesh_importer);
+
+void carve_deleteMesh(struct CarveMeshDescr *mesh_descr);
+
+bool carve_performBooleanOperation(struct CarveMeshDescr *left_mesh,
+                                   struct CarveMeshDescr *right_mesh,
+                                   int operation,
+                                   struct CarveMeshDescr **output_mesh);
+
+void carve_exportMesh(struct CarveMeshDescr *mesh_descr,
+                      CarveMeshExporter *mesh_exporter,
+                      struct ExportMeshData *export_data);
+
+void carve_unionIntersections(struct CarveMeshDescr **left_mesh_r, struct CarveMeshDescr **right_mesh_r);
+
+#ifdef __cplusplus
+}  // extern "C"
+#endif
+
+#endif  // __CARVE_CAPI_H__
diff --git a/extern/carve/carve-util.cc b/extern/carve/carve-util.cc
new file mode 100644
index 00000000000..b1f9ffbb769
--- /dev/null
+++ b/extern/carve/carve-util.cc
@@ -0,0 +1,778 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2014 Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Blender Foundation,
+ *                 Sergey Sharybin
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#include "carve-util.h"
+#include "carve-capi.h"
+
+#include <cfloat>
+#include <carve/csg.hpp>
+#include <carve/csg_triangulator.hpp>
+#include <carve/rtree.hpp>
+
+using carve::csg::Intersections;
+using carve::geom::aabb;
+using carve::geom::RTreeNode;
+using carve::geom3d::Vector;
+using carve::math::Matrix3;
+using carve::mesh::Face;
+using carve::mesh::MeshSet;
+using carve::triangulate::triangulate;
+
+typedef std::map< MeshSet<3>::mesh_t*, RTreeNode<3, Face<3> *> * > RTreeCache;
+typedef std::map< MeshSet<3>::mesh_t*, bool > IntersectCache;
+
+namespace {
+
+// Functions adopted from BLI_math.h to use Carve Vector and Matrix.
+
+void axis_angle_normalized_to_mat3(const Vector &normal,
+                                   const double angle,
+                                   Matrix3 *matrix)
+{
+	double nsi[3], co, si, ico;
+
+	/* now convert this to a 3x3 matrix */
+	co = cos(angle);
+	si = sin(angle);
+
+	ico = (1.0 - co);
+	nsi[0] = normal[0] * si;
+	nsi[1] = normal[1] * si;
+	nsi[2] = normal[2] * si;
+
+	matrix->m[0][0] = ((normal[0] * normal[0]) * ico) + co;
+	matrix->m[0][1] = ((normal[0] * normal[1]) * ico) + nsi[2];
+	matrix->m[0][2] = ((normal[0] * normal[2]) * ico) - nsi[1];
+	matrix->m[1][0] = ((normal[0] * normal[1]) * ico) - nsi[2];
+	matrix->m[1][1] = ((normal[1] * normal[1]) * ico) + co;
+	matrix->m[1][2] = ((normal[1] * normal[2]) * ico) + nsi[0];
+	matrix->m[2][0] = ((normal[0] * normal[2]) * ico) + nsi[1];
+	matrix->m[2][1] = ((normal[1] * normal[2]) * ico) - nsi[0];
+	matrix->m[2][2] = ((normal[2] * normal[2]) * ico) + co;
+}
+
+void axis_angle_to_mat3(const Vector &axis,
+                        const double angle,
+                        Matrix3 *matrix)
+{
+	if (axis.length2() < FLT_EPSILON) {
+		*matrix = Matrix3();
+		return;
+	}
+
+	Vector nor = axis;
+	nor.normalize();
+
+	axis_angle_normalized_to_mat3(nor, angle, matrix);
+}
+
+inline double saacos(double fac)
+{
+	if      (fac <= -1.0) return M_PI;
+	else if (fac >=  1.0) return 0.0;
+	else                  return acos(fac);
+}
+
+bool axis_dominant_v3_to_m3(const Vector &normal,
+                            Matrix3 *matrix)
+{
+	Vector up;
+	Vector axis;
+	double angle;
+
+	up.x = 0.0;
+	up.y = 0.0;
+	up.z = 1.0;
+
+	axis = carve::geom::cross(normal, up);
+	angle = saacos(carve::geom::dot(normal, up));
+
+	if (angle >= FLT_EPSILON) {
+		if (axis.length2() < FLT_EPSILON) {
+			axis[0] = 0.0;
+			axis[1] = 1.0;
+			axis[2] = 0.0;
+		}
+
+		axis_angle_to_mat3(axis, angle, matrix);
+		return true;
+	}
+	else {
+		*matrix = Matrix3();
+		return false;
+	}
+}
+
+void meshset_minmax(const MeshSet<3> *mesh,
+                    Vector *min,
+                    Vector *max)
+{
+	for (uint i = 0; i < mesh->vertex_storage.size(); ++i) {
+		min->x = std::min(min->x, mesh->vertex_storage[i].v.x);
+		min->y = std::min(min->y, mesh->vertex_storage[i].v.y);
+		min->z = std::min(min->z, mesh->vertex_storage[i].v.z);
+		max->x = std::max(max->x, mesh->vertex_storage[i].v.x);
+		max->y = std::max(max->y, mesh->vertex_storage[i].v.y);
+		max->z = std::max(max->z, mesh->vertex_storage[i].v.z);
+	}
+}
+
+}  // namespace
+
+void carve_getRescaleMinMax(const MeshSet<3> *left,
+                            const MeshSet<3> *right,
+                            Vector *min,
+                            Vector *max)
+{
+	min->x = max->x = left->vertex_storage[0].v.x;
+	min->y = max->y = left->vertex_storage[0].v.y;
+	min->z = max->z = left->vertex_storage[0].v.z;
+
+	meshset_minmax(left, min, max);
+	meshset_minmax(right, min, max);
+
+	// Make sure we don't scale object with zero scale.
+	if (std::abs(min->x - max->x) < carve::EPSILON) {
+		min->x = -1.0;
+		max->x = 1.0;
+	}
+	if (std::abs(min->y - max->y) < carve::EPSILON) {
+		min->y = -1.0;
+		max->y = 1.0;
+	}
+	if (std::abs(min->z - max->z) < carve::EPSILON) {
+		min->z = -1.0;
+		max->z = 1.0;
+	}
+}
+
+namespace {
+
+void copyMeshes(const std::vector<MeshSet<3>::mesh_t*> &meshes,
+                std::vector<MeshSet<3>::mesh_t*> *new_meshes)
+{
+	std::vector<MeshSet<3>::mesh_t*>::const_iterator it = meshes.begin();
+	new_meshes->reserve(meshes.size());
+	for (; it != meshes.end(); it++) {
+		MeshSet<3>::mesh_t *mesh = *it;
+		MeshSet<3>::mesh_t *new_mesh = new MeshSet<3>::mesh_t(mesh->faces);
+
+		new_meshes->push_back(new_mesh);
+	}
+}
+
+MeshSet<3> *meshSetFromMeshes(const std::vector<MeshSet<3>::mesh_t*> &meshes)
+{
+	std::vector<MeshSet<3>::mesh_t*> new_meshes;
+
+	copyMeshes(meshes, &new_meshes);
+
+	return new MeshSet<3>(new_meshes);
+}
+
+MeshSet<3> *meshSetFromTwoMeshes(const std::vector<MeshSet<3>::mesh_t*> &left_meshes,
+                                 const std::vector<MeshSet<3>::mesh_t*> &right_meshes)
+{
+	std::vector<MeshSet<3>::mesh_t*> new_meshes;
+
+	copyMeshes(left_meshes, &new_meshes);
+	copyMeshes(right_meshes, &new_meshes);
+
+	return new MeshSet<3>(new_meshes);
+}
+
+bool checkEdgeFaceIntersections_do(Intersections &intersections,
+                                   MeshSet<3>::face_t *face_a,
+                                   MeshSet<3>::edge_t *edge_b)
+{
+	if (intersections.intersects(edge_b, face_a))
+		return true;
+
+	carve::mesh::MeshSet<3>::vertex_t::vector_t _p;
+	if (face_a->simpleLineSegmentIntersection(carve::geom3d::LineSegment(edge_b->v1()->v, edge_b->v2()->v), _p))
+		return true;
+
+	return false;
+}
+
+bool checkEdgeFaceIntersections(Intersections &intersections,
+                                MeshSet<3>::face_t *face_a,
+                                MeshSet<3>::face_t *face_b)
+{
+	MeshSet<3>::edge_t *edge_b;
+
+	edge_b = face_b->edge;
+	do {
+		if (checkEdgeFaceIntersections_do(intersections, face_a, edge_b))
+			return true;
+		edge_b = edge_b->next;
+	} while (edge_b != face_b->edge);
+
+	return false;
+}
+
+inline bool facesAreCoplanar(const MeshSet<3>::face_t *a, const MeshSet<3>::face_t *b)
+{
+	carve::geom3d::Ray temp;
+	// XXX: Find a better definition. This may be a source of problems
+	// if floating point inaccuracies cause an incorrect answer.
+	return !carve::geom3d::planeIntersection(a->plane, b->plane, temp);
+}
+
+bool checkMeshSetInterseciton_do(Intersections &intersections,
+                                 const RTreeNode<3, Face<3> *> *a_node,
+                                 const RTreeNode<3, Face<3> *> *b_node,
+                                 bool descend_a = true)
+{
+	if (!a_node->bbox.intersects(b_node->bbox)) {
+		return false;
+	}
+
+	if (a_node->child && (descend_a || !b_node->child)) {
+		for (RTreeNode<3, Face<3> *> *node = a_node->child; node; node = node->sibling) {
+			if (checkMeshSetInterseciton_do(intersections, node, b_node, false)) {
+				return true;
+			}
+		}
+	}
+	else if (b_node->child) {
+		for (RTreeNode<3, Face<3> *> *node = b_node->child; node; node = node->sibling) {
+			if (checkMeshSetInterseciton_do(intersections, a_node, node, true)) {
+				return true;
+			}
+		}
+	}
+	else {
+		for (size_t i = 0; i < a_node->data.size(); ++i) {
+			MeshSet<3>::face_t *fa = a_node->data[i];
+			aabb<3> aabb_a = fa->getAABB();
+			if (aabb_a.maxAxisSeparation(b_node->bbox) > carve::EPSILON) {
+				continue;
+			}
+
+			for (size_t j = 0; j < b_node->data.size(); ++j) {
+				MeshSet<3>::face_t *fb = b_node->data[j];
+				aabb<3> aabb_b = fb->getAABB();
+				if (aabb_b.maxAxisSeparation(aabb_a) > carve::EPSILON) {
+					continue;
+				}
+
+				std::pair<double, double> a_ra = fa->rangeInDirection(fa->plane.N, fa->edge->vert->v);
+				std::pair<double, double> b_ra = fb->rangeInDirection(fa->plane.N, fa->edge->vert->v);
+				if (carve::rangeSeparation(a_ra, b_ra) > carve::EPSILON) {
+					continue;
+				}
+
+				std::pair<double, double> a_rb = fa->rangeInDirection(fb->plane.N, fb->edge->vert->v);
+				std::pair<double, double> b_rb = fb->rangeInDirection(fb->plane.N, fb->edge->vert->v);
+				if (carve::rangeSeparation(a_rb, b_rb) > carve::EPSILON) {
+					continue;
+				}
+
+				if (!facesAreCoplanar(fa, fb)) {
+					if (checkEdgeFaceIntersections(intersections, fa, fb)) {
+						return true;
+					}
+				}
+			}
+		}
+	}
+
+	return false;
+}
+
+bool checkMeshSetInterseciton(RTreeNode<3, Face<3> *> *rtree_a, RTreeNode<3, Face<3> *> *rtree_b)
+{
+	Intersections intersections;
+	return checkMeshSetInterseciton_do(intersections, rtree_a, rtree_b);
+}
+
+void getIntersectedOperandMeshes(std::vector<MeshSet<3>::mesh_t*> *meshes,
+                                 const MeshSet<3>::aabb_t &otherAABB,
+                                 std::vector<MeshSet<3>::mesh_t*> *operandMeshes,
+                                 RTreeCache *rtree_cache,
+                                 IntersectCache *intersect_cache)
+{
+	std::vector<MeshSet<3>::mesh_t*>::iterator it = meshes->begin();
+	std::vector< RTreeNode<3, Face<3> *> *> meshRTree;
+
+	while (it != meshes->end()) {
+		MeshSet<3>::mesh_t *mesh = *it;
+		bool isAdded = false;
+
+		RTreeNode<3, Face<3> *> *rtree;
+		bool intersects;
+
+		RTreeCache::iterator rtree_found = rtree_cache->find(mesh);
+		if (rtree_found != rtree_cache->end()) {
+			rtree = rtree_found->second;
+		}
+		else {
+			rtree = RTreeNode<3, Face<3> *>::construct_STR(mesh->faces.begin(), mesh->faces.end(), 4, 4);
+			(*rtree_cache)[mesh] = rtree;
+		}
+
+		IntersectCache::iterator intersect_found = intersect_cache->find(mesh);
+		if (intersect_found != intersect_cache->end()) {
+			intersects = intersect_found->second;
+		}
+		else {
+			intersects = rtree->bbox.intersects(otherAABB);
+			(*intersect_cache)[mesh] = intersects;
+		}
+
+		if (intersects) {
+			bool isIntersect = false;
+
+			std::vector<MeshSet<3>::mesh_t*>::iterator operand_it = operandMeshes->begin();
+			std::vector<RTreeNode<3, Face<3> *> *>::iterator tree_it = meshRTree.begin();
+			for (; operand_it!=operandMeshes->end(); operand_it++, tree_it++) {
+				RTreeNode<3, Face<3> *> *operandRTree = *tree_it;
+
+				if (checkMeshSetInterseciton(rtree, operandRTree)) {
+					isIntersect = true;
+					break;
+				}
+			}
+
+			if (!isIntersect) {
+				operandMeshes->push_back(mesh);
+				meshRTree.push_back(rtree);
+
+				it = meshes->erase(it);
+				isAdded = true;
+			}
+		}
+
+		if (!isAdded) {
+			//delete rtree;
+			it++;
+		}
+	}
+
+	std::vector<RTreeNode<3, Face<3> *> *>::iterator tree_it = meshRTree.begin();
+	for (; tree_it != meshRTree.end(); tree_it++) {
+		//delete *tree_it;
+	}
+}
+
+MeshSet<3> *getIntersectedOperand(std::vector<MeshSet<3>::mesh_t*> *meshes,
+                                  const MeshSet<3>::aabb_t &otherAABB,
+                                  RTreeCache *rtree_cache,
+                                  IntersectCache *intersect_cache)
+{
+	std::vector<MeshSet<3>::mesh_t*> operandMeshes;
+	getIntersectedOperandMeshes(meshes, otherAABB, &operandMeshes, rtree_cache, intersect_cache);
+
+	if (operandMeshes.size() == 0)
+		return NULL;
+
+	return meshSetFromMeshes(operandMeshes);
+}
+
+MeshSet<3> *unionIntersectingMeshes(carve::csg::CSG *csg,
+                                    MeshSet<3> *poly,
+                                    const MeshSet<3>::aabb_t &otherAABB)
+{
+	if (poly->meshes.size() <= 1) {
+		return poly;
+	}
+
+	std::vector<MeshSet<3>::mesh_t*> orig_meshes =
+			std::vector<MeshSet<3>::mesh_t*>(poly->meshes.begin(), poly->meshes.end());
+
+	RTreeCache rtree_cache;
+	IntersectCache intersect_cache;
+
+	MeshSet<3> *left = getIntersectedOperand(&orig_meshes,
+	                                         otherAABB,
+	                                         &rtree_cache,
+	                                         &intersect_cache);
+
+	if (!left) {
+		// No maniforlds which intersects another object at all.
+		return poly;
+	}
+
+	while (orig_meshes.size()) {
+		MeshSet<3> *right = getIntersectedOperand(&orig_meshes,
+		                                          otherAABB,
+		                                          &rtree_cache,
+		                                          &intersect_cache);
+
+		if (!right) {
+			// No more intersecting manifolds which intersects other object
+			break;
+		}
+
+		try {
+			if (left->meshes.size()==0) {
+				delete left;
+
+				left = right;
+			}
+			else {
+				MeshSet<3> *result = csg->compute(left, right,
+				                                  carve::csg::CSG::UNION,
+				                                  NULL, carve::csg::CSG::CLASSIFY_EDGE);
+
+				delete left;
+				delete right;
+
+				left = result;
+			}
+		}
+		catch (carve::exception e) {
+			std::cerr << "CSG failed, exception " << e.str() << std::endl;
+
+			MeshSet<3> *result = meshSetFromTwoMeshes(left->meshes, right->meshes);
+
+			delete left;
+			delete right;
+
+			left = result;
+		}
+		catch (...) {
+			delete left;
+			delete right;
+
+			throw "Unknown error in Carve library";
+		}
+	}
+
+	for (RTreeCache::iterator it = rtree_cache.begin();
+	     it != rtree_cache.end();
+	     it++)
+	{
+		delete it->second;
+	}
+
+	// Append all meshes which doesn't have intersection with another operand as-is.
+	if (orig_meshes.size()) {
+		MeshSet<3> *result = meshSetFromTwoMeshes(left->meshes, orig_meshes);
+
+		delete left;
+		left = result;
+	}
+
+	return left;
+}
+
+}  // namespace
+
+// TODO(sergey): This function is to be totally re-implemented to make it
+// more clear what's going on and hopefully optimize it as well.
+void carve_unionIntersections(carve::csg::CSG *csg,
+                              MeshSet<3> **left_r,
+                              MeshSet<3> **right_r)
+{
+	MeshSet<3> *left = *left_r, *right = *right_r;
+
+	if (left->meshes.size() == 1 && right->meshes.size() == 0) {
+		return;
+	}
+
+	MeshSet<3>::aabb_t leftAABB = left->getAABB();
+	MeshSet<3>::aabb_t rightAABB = right->getAABB();;
+
+	left = unionIntersectingMeshes(csg, left, rightAABB);
+	right = unionIntersectingMeshes(csg, right, leftAABB);
+
+	if (left != *left_r) {
+		delete *left_r;
+	}
+
+	if (right != *right_r)
+		delete *right_r;
+
+	*left_r = left;
+	*right_r = right;
+}
+
+static inline void add_newell_cross_v3_v3v3(const Vector &v_prev,
+                                            const Vector &v_curr,
+                                            Vector *n)
+{
+	(*n)[0] += (v_prev[1] - v_curr[1]) * (v_prev[2] + v_curr[2]);
+	(*n)[1] += (v_prev[2] - v_curr[2]) * (v_prev[0] + v_curr[0]);
+	(*n)[2] += (v_prev[0] - v_curr[0]) * (v_prev[1] + v_curr[1]);
+}
+
+// Axis matrix is being set for non-flat ngons only.
+bool carve_checkPolyPlanarAndGetNormal(const std::vector<Vector> &vertices,
+                                       const int verts_per_poly,
+                                       const int *verts_of_poly,
+                                       Matrix3 *axis_matrix_r)
+{
+	if (verts_per_poly == 3) {
+		// Triangles are always planar.
+		return true;
+	}
+	else if (verts_per_poly == 4) {
+		// Presumably faster than using generig n-gon check for quads.
+
+		const Vector &v1 = vertices[verts_of_poly[0]],
+		             &v2 = vertices[verts_of_poly[1]],
+		             &v3 = vertices[verts_of_poly[2]],
+		             &v4 = vertices[verts_of_poly[3]];
+
+		Vector vec1, vec2, vec3, cross;
+
+		vec1 = v2 - v1;
+		vec2 = v4 - v1;
+		vec3 = v3 - v1;
+
+		cross = carve::geom::cross(vec1, vec2);
+
+		double production = carve::geom::dot(cross, vec3);
+
+		// TODO(sergey): Check on whether we could have length-independent
+		// magnitude here.
+		double magnitude = 1e-3 * cross.length2();
+
+		return fabs(production) < magnitude;
+	}
+	else {
+		const Vector *vert_prev = &vertices[verts_of_poly[verts_per_poly - 1]];
+		const Vector *vert_curr = &vertices[verts_of_poly[0]];
+
+		Vector normal = carve::geom::VECTOR(0.0, 0.0, 0.0);
+		for (int i = 0; i < verts_per_poly; i++) {
+			add_newell_cross_v3_v3v3(*vert_prev, *vert_curr, &normal);
+			vert_prev = vert_curr;
+			vert_curr = &vertices[verts_of_poly[(i + 1) % verts_per_poly]];
+		}
+
+		if (normal.length2() < FLT_EPSILON) {
+			// Degenerated face, couldn't triangulate properly anyway.
+			return true;
+		}
+		else {
+			double magnitude = normal.length2();
+
+			normal.normalize();
+			axis_dominant_v3_to_m3(normal, axis_matrix_r);
+
+			Vector first_projected = *axis_matrix_r * vertices[verts_of_poly[0]];
+			double min_z = first_projected[2], max_z = first_projected[2];
+
+			for (int i = 1; i < verts_per_poly; i++) {
+				const Vector &vertex = vertices[verts_of_poly[i]];
+				Vector projected = *axis_matrix_r * vertex;
+				if (projected[2] < min_z) {
+					min_z = projected[2];
+				}
+				if (projected[2] > max_z) {
+					max_z = projected[2];
+				}
+			}
+
+			if (std::abs(min_z - max_z) > FLT_EPSILON * magnitude) {
+				return false;
+			}
+		}
+
+		return true;
+	}
+
+	return false;
+}
+
+namespace {
+
+int triangulateNGon_carveTriangulator(const std::vector<Vector> &vertices,
+                                      const int verts_per_poly,
+                                      const int *verts_of_poly,
+                                      const Matrix3 &axis_matrix,
+                                      std::vector<carve::triangulate::tri_idx> *triangles)
+{
+	// Project vertices to 2D plane.
+	Vector projected;
+	std::vector<carve::geom::vector<2> > poly_2d;
+	poly_2d.reserve(verts_per_poly);
+	for (int i = 0; i < verts_per_poly; ++i) {
+		projected = axis_matrix * vertices[verts_of_poly[i]];
+		poly_2d.push_back(carve::geom::VECTOR(projected[0], projected[1]));
+	}
+
+	carve::triangulate::triangulate(poly_2d, *triangles);
+
+	return triangles->size();
+}
+
+int triangulateNGon_importerTriangulator(struct ImportMeshData *import_data,
+                                         CarveMeshImporter *mesh_importer,
+                                         const std::vector<Vector> &vertices,
+                                         const int verts_per_poly,
+                                         const int *verts_of_poly,
+                                         const Matrix3 &axis_matrix,
+                                         std::vector<carve::triangulate::tri_idx> *triangles)
+{
+	typedef float Vector2D[2];
+	typedef unsigned int Triangle[3];
+
+	// Project vertices to 2D plane.
+	Vector2D *poly_2d = new Vector2D[verts_per_poly];
+	Vector projected;
+	for (int i = 0; i < verts_per_poly; ++i) {
+		projected = axis_matrix * vertices[verts_of_poly[i]];
+		poly_2d[i][0] = projected[0];
+		poly_2d[i][1] = projected[1];
+	}
+
+	Triangle *api_triangles = new Triangle[verts_per_poly - 2];
+	int num_triangles =
+		mesh_importer->triangulate2DPoly(import_data,
+		                                 poly_2d,
+		                                 verts_per_poly,
+		                                 api_triangles);
+
+	triangles->reserve(num_triangles);
+	for (int i = 0; i < num_triangles; ++i) {
+		triangles->push_back(
+			carve::triangulate::tri_idx(api_triangles[i][0],
+			                            api_triangles[i][1],
+			                            api_triangles[i][2]));
+	}
+
+	delete poly_2d;
+	delete api_triangles;
+
+	return num_triangles;
+}
+
+}  // namespace
+
+int carve_triangulatePoly(struct ImportMeshData *import_data,
+                          CarveMeshImporter *mesh_importer,
+                          int poly_index,
+                          int start_loop_index,
+                          const std::vector<Vector> &vertices,
+                          const int verts_per_poly,
+                          const int *verts_of_poly,
+                          const Matrix3 &axis_matrix,
+                          std::vector<int> *face_indices,
+                          std::vector<int> *orig_loop_index_map,
+                          std::vector<int> *orig_poly_index_map)
+{
+	int num_triangles = 0;
+
+	assert(verts_per_poly > 3);
+
+	if (verts_per_poly == 4) {
+		// Quads we triangulate by 1-3 diagonal, it is an original behavior
+		// of boolean modifier.
+		//
+		// TODO(sergey): Consider using shortest diagonal here. However
+		// display code in Blende use static 1-3 split, so some experiments
+		// are needed here.
+		face_indices->push_back(3);
+		face_indices->push_back(verts_of_poly[0]);
+		face_indices->push_back(verts_of_poly[1]);
+		face_indices->push_back(verts_of_poly[2]);
+
+		orig_loop_index_map->push_back(start_loop_index);
+		orig_loop_index_map->push_back(start_loop_index + 1);
+		orig_loop_index_map->push_back(-1);
+		orig_poly_index_map->push_back(poly_index);
+
+		face_indices->push_back(3);
+		face_indices->push_back(verts_of_poly[0]);
+		face_indices->push_back(verts_of_poly[2]);
+		face_indices->push_back(verts_of_poly[3]);
+
+		orig_loop_index_map->push_back(-1);
+		orig_loop_index_map->push_back(start_loop_index + 2);
+		orig_loop_index_map->push_back(start_loop_index + 3);
+		orig_poly_index_map->push_back(poly_index);
+
+		num_triangles = 2;
+	}
+	else {
+		std::vector<carve::triangulate::tri_idx> triangles;
+		triangles.reserve(verts_per_poly - 2);
+
+		// Make triangulator callback optional so we could do some tests
+		// in the future.
+		if (mesh_importer->triangulate2DPoly) {
+			num_triangles =
+				triangulateNGon_importerTriangulator(import_data,
+				                                     mesh_importer,
+				                                     vertices,
+				                                     verts_per_poly,
+				                                     verts_of_poly,
+				                                     axis_matrix,
+				                                     &triangles);
+		}
+		else {
+			num_triangles =
+				triangulateNGon_carveTriangulator(vertices,
+				                                  verts_per_poly,
+				                                  verts_of_poly,
+				                                  axis_matrix,
+				                                  &triangles);
+		}
+
+		for (int i = 0; i < triangles.size(); ++i) {
+			int v1 = triangles[i].c,
+			    v2 = triangles[i].b,
+			    v3 = triangles[i].a;
+
+			// Sanity check of the triangle.
+			assert(v1 != v2);
+			assert(v1 != v3);
+			assert(v2 != v3);
+			assert(v1 < verts_per_poly);
+			assert(v2 < verts_per_poly);
+			assert(v3 < verts_per_poly);
+
+			face_indices->push_back(3);
+			face_indices->push_back(verts_of_poly[v3]);
+			face_indices->push_back(verts_of_poly[v2]);
+			face_indices->push_back(verts_of_poly[v1]);
+
+#define CHECK_TRIANGLE_LOOP_INDEX(v1, v2) \
+	{ \
+		if (v2 == v1 + 1) { \
+			orig_loop_index_map->push_back(start_loop_index + v1); \
+		} \
+		else { \
+			orig_loop_index_map->push_back(-1); \
+		} \
+	} (void) 0
+
+			CHECK_TRIANGLE_LOOP_INDEX(v1, v2);
+			CHECK_TRIANGLE_LOOP_INDEX(v2, v3);
+			CHECK_TRIANGLE_LOOP_INDEX(v3, v1);
+
+#undef CHECK_TRIANGLE_LOOP_INDEX
+
+			orig_poly_index_map->push_back(poly_index);
+		}
+	}
+
+	return num_triangles;
+}
diff --git a/extern/carve/carve-util.h b/extern/carve/carve-util.h
new file mode 100644
index 00000000000..07743de2d16
--- /dev/null
+++ b/extern/carve/carve-util.h
@@ -0,0 +1,122 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2014 Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Blender Foundation,
+ *                 Sergey Sharybin
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#ifndef __CARVE_UTIL_H__
+#define __CARVE_UTIL_H__
+
+#include <carve/csg.hpp>
+#include <carve/geom3d.hpp>
+#include <carve/interpolator.hpp>
+#include <carve/mesh.hpp>
+
+#include "carve-capi.h"
+
+void carve_getRescaleMinMax(const carve::mesh::MeshSet<3> *left,
+                            const carve::mesh::MeshSet<3> *right,
+                            carve::geom3d::Vector *min,
+                            carve::geom3d::Vector *max);
+
+void carve_unionIntersections(carve::csg::CSG *csg,
+                              carve::mesh::MeshSet<3> **left_r,
+                              carve::mesh::MeshSet<3> **right_r);
+
+bool carve_checkPolyPlanarAndGetNormal(const std::vector<carve::geom3d::Vector> &vertices,
+                                       const int verts_per_poly,
+                                       const int *verts_of_poly,
+                                       carve::math::Matrix3 *axis_matrix_r);
+
+int carve_triangulatePoly(struct ImportMeshData *import_data,
+                          CarveMeshImporter *mesh_importer,
+                          int poly_index,
+                          int start_loop_index,
+                          const std::vector<carve::geom3d::Vector> &vertices,
+                          const int verts_per_poly,
+                          const int *verts_of_poly,
+                          const carve::math::Matrix3 &axis_matrix,
+                          std::vector<int> *face_indices,
+                          std::vector<int> *orig_loop_index_map,
+                          std::vector<int> *orig_poly_index_map);
+
+namespace carve {
+	namespace interpolate {
+
+		template<typename attr_t>
+		class VertexAttr : public Interpolator {
+		public:
+			typedef const meshset_t::vertex_t *key_t;
+
+		protected:
+			typedef std::unordered_map<key_t, attr_t> attrmap_t;
+
+			attrmap_t attrs;
+
+			virtual void resultFace(const carve::csg::CSG &csg,
+			                        const meshset_t::face_t *new_face,
+			                        const meshset_t::face_t *orig_face,
+			                        bool flipped)
+			{
+				typedef meshset_t::face_t::const_edge_iter_t const_edge_iter_t;
+				for (const_edge_iter_t new_edge_iter = new_face->begin();
+					 new_edge_iter != new_face->end();
+					 ++new_edge_iter)
+				{
+					typename attrmap_t::const_iterator found =
+						attrs.find(new_edge_iter->vert);
+					if (found == attrs.end()) {
+						for (const_edge_iter_t orig_edge_iter = orig_face->begin();
+							 orig_edge_iter != orig_face->end();
+							 ++orig_edge_iter)
+						{
+							if ((orig_edge_iter->vert->v - new_edge_iter->vert->v).length2() < 1e-5) {
+								attrs[new_edge_iter->vert] = attrs[orig_edge_iter->vert];
+							}
+						}
+					}
+				}
+			}
+
+		public:
+			bool hasAttribute(const meshset_t::vertex_t *v) {
+				return attrs.find(v) != attrs.end();
+			}
+
+			const attr_t &getAttribute(const meshset_t::vertex_t *v, const attr_t &def = attr_t()) {
+				typename attrmap_t::const_iterator found = attrs.find(v);
+				if (found != attrs.end()) {
+					return found->second;
+				}
+				return def;
+			}
+
+			void setAttribute(const meshset_t::vertex_t *v, const attr_t &attr) {
+				attrs[v] = attr;
+			}
+		};
+
+	}  // namespace interpolate
+}  // namespace carve
+
+#endif  // __CARVE_UTIL_H__
diff --git a/extern/carve/include/carve/interpolator.hpp b/extern/carve/include/carve/interpolator.hpp
index f8162e52f39..b88a90dd18d 100644
--- a/extern/carve/include/carve/interpolator.hpp
+++ b/extern/carve/include/carve/interpolator.hpp
@@ -219,7 +219,7 @@ namespace carve {
           interpolator->edgeDivision(csg, orig_edge, orig_edge_idx, v1, v2);
         }
 
-        Hook(Interpolator *_interpolator, const carve::csg::CSG &_csg) : interpolator(_interpolator), csg(_csg) {
+        Hook(Interpolator *_interpolator, const carve::csg::CSG &_csg) : csg(_csg), interpolator(_interpolator) {
         }
 
         virtual ~Hook() {
diff --git a/extern/carve/patches/interpolator_reorder.patch b/extern/carve/patches/interpolator_reorder.patch
new file mode 100644
index 00000000000..867297fef7d
--- /dev/null
+++ b/extern/carve/patches/interpolator_reorder.patch
@@ -0,0 +1,12 @@
+diff -r e82d852e4fb0 include/carve/interpolator.hpp
+--- a/include/carve/interpolator.hpp	Wed Jan 15 13:16:14 2014 +1100
++++ b/include/carve/interpolator.hpp	Fri Jan 31 18:55:05 2014 +0600
+@@ -219,7 +219,7 @@
+           interpolator->edgeDivision(csg, orig_edge, orig_edge_idx, v1, v2);
+         }
+ 
+-        Hook(Interpolator *_interpolator, const carve::csg::CSG &_csg) : interpolator(_interpolator), csg(_csg) {
++        Hook(Interpolator *_interpolator, const carve::csg::CSG &_csg) : csg(_csg), interpolator(_interpolator) {
+         }
+ 
+         virtual ~Hook() {
diff --git a/extern/carve/patches/series b/extern/carve/patches/series
index d579cdaf277..30937b4b9cf 100644
--- a/extern/carve/patches/series
+++ b/extern/carve/patches/series
@@ -5,3 +5,4 @@ mingw.patch
 gcc46.patch
 clang_is_heap_fix.patch
 strict_flags.patch
+interpolator_reorder.patch