std::vector<DeletedFacet> InterfaceGeometry::BuildIntersection(size_t *nbCreated) {
	mApp->changedSinceSave = true;
	std::vector<InterfaceVertex> newVertices;
	std::vector<IntersectFacet> selectedFacets;
	std::vector<DeletedFacet> deletedFacetList;
	*nbCreated = 0; //Total number of new facets created

	//Populate selected facets
	for (size_t i = 0; i < sh.nbFacet; i++) {
		if (facets[i]->selected) {
			IntersectFacet facet;
			facet.id = i;
			facet.f = facets[i];
			facet.intersectionPointId.resize(facet.f->sh.nbIndex);
			selectedFacets.push_back(facet);
		}
	}
	for (size_t i = 0; i < selectedFacets.size(); i++) {
		InterfaceFacet* f1 = selectedFacets[i].f;
		for (size_t j = 0; j < selectedFacets.size(); j++) {
			InterfaceFacet* f2 = selectedFacets[j].f;
			if (i != j) {
				size_t c1, c2, l;
				if (!GetCommonEdges(f1, f2, &c1, &c2, &l)) {
					for (size_t index = 0; index < f1->sh.nbIndex; index++) { //Go through all indexes of edge-finding facet
						InterfaceVertex intersectionPoint;
						InterfaceVertex base = vertices3[f1->indices[index]];
						Vector3d side = vertices3[f1->GetIndex(index + 1)] - base;

						if (IntersectingPlaneWithLine(base, side, f2->sh.O, f2->sh.N, &intersectionPoint, true)) {
							Vector2d projected = ProjectVertex(intersectionPoint, f2->sh.U, f2->sh.V, f2->sh.O);
							bool inPoly = IsInPoly(projected.u, projected.v, f2->vertices2);
							bool onEdge = IsOnPolyEdge(projected.u, projected.v, f2->vertices2, 1E-6);
							//onEdge = false;
							if (inPoly || onEdge) {
								//Intersection found. First check if we already created this point
								int foundId = -1;
								for (size_t v = 0; foundId == -1 && v < newVertices.size(); v++) {
									if (IsZero((newVertices[v] - intersectionPoint).Norme()))
										foundId = (int)v;
								}
								IntersectPoint newPoint{}, newPointOtherFacet{};
								if (foundId == -1) { //Register new intersection point
									newPoint.vertexId = newPointOtherFacet.vertexId = sh.nbVertex + newVertices.size();

									intersectionPoint.selected = false;
									newVertices.push_back(intersectionPoint);
								}
								else { //Refer to existing intersection point
									newPoint.vertexId = newPointOtherFacet.vertexId = foundId + sh.nbVertex;
								}
								newPoint.withFacetId = j;
								selectedFacets[i].intersectionPointId[index].push_back(newPoint);
								newPointOtherFacet.withFacetId = i; //With my edge
								if (!onEdge) selectedFacets[j].intersectingPoints.push_back(newPointOtherFacet); //Other facet's plane intersected
							}
						}
					}
				}
			}
		}
	}

	vertices3.insert(vertices3.end(), newVertices.begin(), newVertices.end());
	sh.nbVertex += newVertices.size();

	UnselectAll();
	for (auto & selectedFacet : selectedFacets) {
		std::vector<std::vector<EdgePoint>> clipPaths;
		InterfaceFacet *f = selectedFacet.f;
		for (size_t vertexId = 0; vertexId < selectedFacet.f->sh.nbIndex; vertexId++) { //Go through indices
			//testPath.push_back(f->indices[vertexId]);
			for (size_t v = 0; v < selectedFacet.intersectionPointId[vertexId].size(); v++) { //If there are intersection points on this edge, go through them
				//Check if not the end of an already registered clipping path
				bool found = false;
				for (size_t i = 0; found == false && i < clipPaths.size(); i++) {
					found = clipPaths[i].back().vertexId == selectedFacet.intersectionPointId[vertexId][v].vertexId;
				}
				if (!found) { //Register a new clip path
					std::vector<EdgePoint> path;
					EdgePoint p{};
					p.vertexId = selectedFacet.intersectionPointId[vertexId][v].vertexId;
					p.onEdge = (int)vertexId;
					path.push_back(p); //Register intersection point
					size_t searchId = selectedFacet.intersectionPointId[vertexId][v].vertexId; //Current point, by global Id
					size_t searchFacetId = selectedFacet.intersectionPointId[vertexId][v].withFacetId; // Current facet with which we intersected
					//v++;
					int foundId, foundId2;
					do {
						foundId = -1;
						for (size_t p1 = 0; foundId == -1 && p1 < selectedFacet.intersectingPoints.size(); p1++) { //Get the next intersection point with same facet
							if (searchId == p1) continue;
							foundId = ((selectedFacet.intersectingPoints[p1].withFacetId == searchFacetId) && (selectedFacet.intersectingPoints[p1].vertexId != searchId)) ? (int)p1 : -1;
						}
						if (foundId != -1) {
							EdgePoint point{};
                            point.vertexId = selectedFacet.intersectingPoints[foundId].vertexId;
                            point.onEdge = -1;
							path.push_back(point);
							//Get next point which is the same
							searchId = selectedFacet.intersectingPoints[foundId].vertexId;
							searchFacetId = selectedFacet.intersectingPoints[foundId].withFacetId;
							foundId2 = -1;
							for (size_t p2 = 0; foundId2 == -1 && p2 < selectedFacet.intersectingPoints.size(); p2++) { //Search next intersection point which is same vertex
								if (p2 == foundId) continue;
								foundId2 = ((selectedFacet.intersectingPoints[p2].vertexId == searchId) && (selectedFacet.intersectingPoints[p2].withFacetId != searchFacetId)) ? (int)p2 : -1;
							}
							if (foundId2 != -1) {
								searchFacetId = selectedFacet.intersectingPoints[foundId2].withFacetId;
								searchId = foundId2;
							}
						}
					} while (foundId != -1 && foundId2 != -1);
					//No more intersection points on the middle of the facet. Need to find closing point, which is on an edge
					foundId = -1;
					for (size_t v2 = v; foundId == -1 && v2 < selectedFacet.intersectionPointId[vertexId].size(); v2++) { //Check if on same edge
						if (selectedFacet.intersectionPointId[vertexId][v2].withFacetId == searchFacetId && selectedFacet.intersectionPointId[vertexId][v2].vertexId != path.front().vertexId) foundId = (int)selectedFacet.intersectionPointId[vertexId][v2].vertexId;
					}
					if (foundId != -1) {
						EdgePoint point{};
                        point.vertexId = foundId;
                        point.onEdge = (int)vertexId;
						path.push_back(point); //Found on same edge, close
					}
					else { //Search on other edges
						for (size_t v3 = 0; foundId == -1 && v3 < selectedFacet.f->sh.nbIndex; v3++) {
							if (v3 == vertexId) continue; //Already checked on same edge
							for (size_t v2 = v; foundId == -1 && v2 < selectedFacet.intersectionPointId[v3].size(); v2++) {
								if (selectedFacet.intersectionPointId[v3][v2].withFacetId == searchFacetId) {
									foundId = (int)selectedFacet.intersectionPointId[v3][v2].vertexId;
									EdgePoint point{};
                                    point.vertexId = foundId;
                                    point.onEdge = (int)v3;
									path.push_back(point);
								}
							}
						}
					}
					if (path.size() > 1) clipPaths.push_back(path);
				}
			}
		}

		if (!clipPaths.empty()) {
			//Construct clipped facet, having a selected vertex
			std::vector<bool> isIndexSelected(f->sh.nbIndex);
			for (size_t v = 0; v < f->sh.nbIndex; v++) {
				isIndexSelected[v] = vertices3[f->indices[v]].selected; //Make a copy, we don't want to deselect vertices
			}
			size_t nbNewfacet = 0; //Number of new facets created for this particular clipping path
			for (size_t v = 0; v < f->sh.nbIndex; v++) {
				size_t currentVertex = v;
				if (isIndexSelected[currentVertex]) {
					//Restore visited state for all clip paths
					for (auto & clipPath : clipPaths) {
						for (auto & j : clipPath) {
							j.visited = false;
						}
					}
					std::vector<size_t> clipPath;
					nbNewfacet++;
					do { //Build points of facet
						clipPath.push_back(f->indices[currentVertex]);
						isIndexSelected[currentVertex] = false;
						//Get closest path end
						double minDist = 9E99;
						int clipId = -1;
						bool front;

						for (size_t i = 0; i < clipPaths.size(); i++) {
							if (clipPaths[i].front().onEdge == currentVertex && clipPaths[i].back().onEdge != -1) { //If a full clipping path is found on the scanned edge, go through it
								double d = (vertices3[f->indices[currentVertex]] - vertices3[clipPaths[i].front().vertexId]).Norme();
								if (d < minDist) {
									minDist = d;
									clipId = (int)i;
									front = true;
								}
							}
							if (clipPaths[i].back().onEdge == currentVertex && clipPaths[i].front().onEdge != -1) { //If a full clipping path is found on the scanned edge, go through it
								double d = (vertices3[f->indices[currentVertex]] - vertices3[clipPaths[i].back().vertexId]).Norme();
								if (d < minDist) {
									minDist = d;
									clipId = (int)i;
									front = false;
								}
							}
						}
						if (clipId != -1) {
							if ((front && !clipPaths[clipId].front().visited) || (!front && !clipPaths[clipId].back().visited)) {
								for (int cp = front ? 0 : (int)clipPaths[clipId].size() - 1; cp >= 0 && cp < (int)clipPaths[clipId].size(); cp += front ? 1 : -1) {
									clipPath.push_back(clipPaths[clipId][cp].vertexId);
									clipPaths[clipId][cp].visited = true;
								}
								currentVertex = front ? clipPaths[clipId].back().onEdge : clipPaths[clipId].front().onEdge;
							}
						}
						currentVertex = (currentVertex + 1) % f->sh.nbIndex;
					} while (currentVertex != v);
					if (clipPath.size() > 2) {
						auto *fac = new InterfaceFacet(clipPath.size());
						for (size_t i = 0; i < clipPath.size(); i++)
                            fac->indices[i] = clipPath[i];
                        fac->selected = true;

						if (nbNewfacet == 1) {//replace original
							DeletedFacet df{};
							df.f = facets[selectedFacet.id];
							df.ori_pos = selectedFacet.id;
							df.replaceOri = true;
							deletedFacetList.push_back(df);
							facets[selectedFacet.id] = fac; //replace original
						}
						else {//create new
                            facets.push_back(fac);
                            sh.nbFacet++;
							(*nbCreated)++;
						}
					}
				}
			}
		}

	Rebuild();
	return deletedFacetList;
}