Commit 8259ea87 by Pierre Kraemer

update several Algo/Geometry functions

parent f0e1b718
 ... ... @@ -25,7 +25,6 @@ #ifndef __ALGO_GEOMETRY_BOUNDINGBOX_H__ #define __ALGO_GEOMETRY_BOUNDINGBOX_H__ #include "Geometry/basic.h" #include "Geometry/bounding_box.h" #include "Topology/generic/attributeHandler.h" #include "Topology/generic/traversor/traversorCell.h" ... ... @@ -43,9 +42,7 @@ template Geom::BoundingBox computeBoundingBox(typename PFP::MAP& map, const VertexAttribute& position) { Geom::BoundingBox bb ; TraversorV t(map) ; for(Dart d = t.begin(); d != t.end(); d = t.next()) bb.addPoint(position[d]) ; foreach_cell(map, [&] (Vertex v) { bb.addPoint(position[v]) ; }); return bb ; } ... ...
 ... ... @@ -41,7 +41,7 @@ namespace Geometry * Test if an edge bounded by 2 faces is convex or concave */ template bool isEdgeConvexe(typename PFP::MAP& map, Dart d, const VertexAttribute& position) ; bool isEdgeConvex(typename PFP::MAP& map, Edge e, const VertexAttribute& position) ; } // namespace Geometry ... ...
 ... ... @@ -40,14 +40,14 @@ namespace Geometry { template bool isEdgeConvexe(typename PFP::MAP& map, Dart d, const VertexAttribute& position) bool isEdgeConvex(typename PFP::MAP& map, Edge e, const VertexAttribute& position) { typedef typename PFP::VEC3 VEC3 ; const VEC3 n = faceNormal(map, d, position); const VEC3 e = vectorOutOfDart(map, map.phi1(map.phi2(d)), position) ; const VEC3 n = faceNormal(map, e.dart, position); const VEC3 ee = vectorOutOfDart(map, map.phi1(map.phi2(e.dart)), position) ; if((e * n) > 0) if((ee * n) < 0) return true; else return false; ... ...
 ... ... @@ -37,38 +37,36 @@ namespace Geometry /** * compute squared distance from point to the plane of a planar face * @param map the map * @param d a dart of the face * @param f a face * @param position the vertex attribute storing positions * @param P the point * @return the squared distance to tha plane * @return the squared distance to the plane */ template typename PFP::REAL squaredDistancePoint2FacePlane(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const VEC3& P) ; typename PFP::REAL squaredDistancePoint2FacePlane(typename PFP::MAP& map, Face f, const VertexAttribute& position, const VEC3& P) ; /** * compute squared distance from point to face (assuming face is convex) * Algo: min distance of each subtriangle of face (not optimum ?) * @param map the map * @param d a dart of the face * @param f a face * @param position the vertex attribute storing positions * @param P the point * @return the squared distance */ template typename PFP::REAL squaredDistancePoint2Face(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const VEC3& P) ; typename PFP::REAL squaredDistancePoint2Face(typename PFP::MAP& map, Face f, const VertexAttribute& position, const VEC3& P) ; /** * compute squared distance from point to an edge * @param map the map * @param d a dart of the edge * @param e an edge * @param position the vertex attribute storing positions * @param P the point * @return the squared distance */ template typename PFP::REAL squaredDistancePoint2Edge(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const VEC3& P) ; template bool isPlanar(typename PFP::MAP& map, Dart d, const VertexAttribute& position); typename PFP::REAL squaredDistancePoint2Edge(typename PFP::MAP& map, Edge e, const VertexAttribute& position, const VEC3& P) ; } // namespace Geometry ... ...
 ... ... @@ -23,8 +23,6 @@ *******************************************************************************/ #include "Geometry/distances.h" #include "Algo/Geometry/normal.h" #include namespace CGoGN { ... ... @@ -36,64 +34,50 @@ namespace Geometry { template bool isPlanar(typename PFP::MAP& map, Dart d, const VertexAttribute& position) typename PFP::REAL squaredDistancePoint2FacePlane(typename PFP::MAP& map, Face f, const VertexAttribute& position, const VEC3& P) { if (map.phi<111>(d)==d) return true; Vertex v(f.dart); const typename PFP::VEC3& A = position[v]; v.dart = map.phi1(v.dart); const typename PFP::VEC3& B = position[v]; v.dart = map.phi1(v.dart); const typename PFP::VEC3& C = position[v]; typename PFP::VEC3 No = Algo::Surface::Geometry::triangleNormal(map, d, position) ; Dart e = map.phi<11>(d); while (e != d) { typename PFP::VEC3 n = Algo::Surface::Geometry::triangleNormal(map, e, position) ; e = map.phi1(e); if (e!=d) e = map.phi1(e); } return Geom::squaredDistancePoint2TrianglePlane(P, A, B, C); } template typename PFP::REAL squaredDistancePoint2FacePlane(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const VEC3& P) { const typename PFP::VEC3& A = position[d]; d = map.phi1(d); const typename PFP::VEC3& B = position[d]; d = map.phi1(d); const typename PFP::VEC3& C = position[d]; return Geom::squaredDistancePoint2TrianglePlane(P,A,B,C); } template typename PFP::REAL squaredDistancePoint2Face(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const VEC3& P) typename PFP::REAL squaredDistancePoint2Face(typename PFP::MAP& map, Face f, const VertexAttribute& position, const VEC3& P) { typedef typename PFP::REAL REAL; const typename PFP::VEC3& A = position[d]; const typename PFP::VEC3& A = position[f.dart]; Dart d = map.phi1(f.dart); Dart e = map.phi1(d); Dart f = map.phi1(e); REAL dist2 = Geom::squaredDistancePoint2Triangle(P,A,position[e],position[f]); e = f; f = map.phi1(e); REAL dist2 = Geom::squaredDistancePoint2Triangle(P, A, position[d], position[e]); d = e; e = map.phi1(d); while (f != d) while (e != f.dart) { REAL d2 = Geom::squaredDistancePoint2Triangle(P,A,position[e],position[f]); REAL d2 = Geom::squaredDistancePoint2Triangle(P, A, position[d], position[e]); if (d2 < dist2) dist2 = d2; e = f; f = map.phi1(e); d = e; e = map.phi1(d); } return dist2; } template typename PFP::REAL squaredDistancePoint2Edge(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const VEC3& P) typename PFP::REAL squaredDistancePoint2Edge(typename PFP::MAP& map, Edge e, const VertexAttribute& position, const VEC3& P) { const typename PFP::VEC3& A = position[d]; typename PFP::VEC3& AB = position[map.phi1(d)]-A; typename PFP::REAL AB2 = AB*AB; return Geom::squaredDistanceSeg2Point(A,AB,AB2,P) ; const typename PFP::VEC3& A = position[e.dart]; typename PFP::VEC3 AB = position[map.phi1(e.dart)] - A; typename PFP::REAL AB2 = AB * AB; return Geom::squaredDistanceSeg2Point(A, AB, AB2, P) ; } } // namespace Geometry ... ...
 ... ... @@ -24,7 +24,7 @@ #include "Geometry/basic.h" #include "Algo/Geometry/normal.h" #include "Topology/generic/traversorCell.h" #include "Topology/generic/traversor/traversorCell.h" namespace CGoGN { ... ...
 ... ... @@ -45,61 +45,61 @@ namespace Geometry /** * test if the volume is convex * @param the map * @param a dart of the volume * @param a volume * @param true if the faces of the volume must be in CCW order (default=true) */ template bool isConvex(typename PFP::MAP& map, Dart d, const VertexAttribute& positions, bool CCW, unsigned int thread = 0); bool isConvex(typename PFP::MAP& map, Vol v, const VertexAttribute& positions, bool CCW, unsigned int thread = 0); /** * test if a point is inside a volume * @param map the map * @param d a dart defining the volume * @param a volume * @param the point */ template bool isPointInVolume(typename PFP::MAP& map, Dart d, const VertexAttribute& positions, const typename PFP::VEC3& point); bool isPointInVolume(typename PFP::MAP& map, Vol v, const VertexAttribute& positions, const typename PFP::VEC3& point); /** * test if a point is inside a volume * @param map the map * @param d a dart defining a convex volume * @param a convex volume * @param the point */ template bool isPointInConvexVolume(typename PFP::MAP& map, Dart d, const VertexAttribute& positions, const typename PFP::VEC3& point, bool CCW = true); bool isPointInConvexVolume(typename PFP::MAP& map, Vol v, const VertexAttribute& positions, const typename PFP::VEC3& point, bool CCW = true); /** * test if a point is inside a face in a plane * test if a point is inside a face * @param map the map * @param d a dart defining the face * @param a face * @param the point */ template bool isPointInConvexFace2D(typename PFP::MAP& map, Dart d, const VertexAttribute& positions, const typename PFP::VEC3& point, bool CCW = true); bool isPointInConvexFace(typename PFP::MAP& map, Face f, const VertexAttribute& positions, const typename PFP::VEC3& point, bool CCW); /** * test if a point is inside a face * test if a point is inside a face in a plane * @param map the map * @param d a dart defining the face * @param a face * @param the point */ template bool isPointInConvexFace(typename PFP::MAP& map, Dart d, const VertexAttribute& positions, const typename PFP::VEC3& point, bool CCW); bool isPointInConvexFace2D(typename PFP::MAP& map, Face f, const VertexAttribute& positions, const typename PFP::VEC3& point, bool CCW = true); /** * test if a point is on an edge * @param map the map * @param d a dart defining the edge * @param an edge * @param the point */ template bool isPointOnEdge(typename PFP::MAP& map, Dart d, const VertexAttribute& positions, const typename PFP::VEC3& point); bool isPointOnEdge(typename PFP::MAP& map, Edge e, const VertexAttribute& positions, const typename PFP::VEC3& point); /** * test if a point is on an half-edge defined by a dart * @param map the map * @param d a dart defining the edge * @param a Dart * @param the point */ template ... ... @@ -108,23 +108,22 @@ bool isPointOnHalfEdge(typename PFP::MAP& map, Dart d, const VertexAttribute bool isPointOnVertex(typename PFP::MAP& map, Dart d, const VertexAttribute& positions, const typename PFP::VEC3& point); bool isPointOnVertex(typename PFP::MAP& map, Vertex v, const VertexAttribute& positions, const typename PFP::VEC3& point); /** * test if a face is intersecting or totally included in a tetrahedron * TODO to test * @param map the map * @param d a dart defining the face * @param f a face * @param the points of the tetra (0,1,2) the first face and (3) the last point of the tetra, in well oriented order * @param true if the faces of the tetra are in CCW order (default=true) */ template bool isConvexFaceInOrIntersectingTetrahedron(typename PFP::MAP& map, Dart d, const VertexAttribute& positions, const typename PFP::VEC3 points[4], bool CCW); bool isConvexFaceInOrIntersectingTetrahedron(typename PFP::MAP& map, Face f, const VertexAttribute& positions, const typename PFP::VEC3 points[4], bool CCW); } // namespace Geometry ... ...
 ... ... @@ -66,7 +66,7 @@ bool isConvex(typename PFP::MAP& map, Vol v, const VertexAttribute bool isPointInVolume(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const typename PFP::VEC3& point) bool isPointInVolume(typename PFP::MAP& map, Vol v, const VertexAttribute& position, const typename PFP::VEC3& point) { typedef typename PFP::VEC3 VEC3; ... ... @@ -75,74 +75,79 @@ bool isPointInVolume(typename PFP::MAP& map, Dart d, const VertexAttribute interPrec; interPrec.reserve(16); std::vector visitedFaces; // Faces that are traversed std::vector visitedFaces; // Faces that are traversed visitedFaces.reserve(64); visitedFaces.push_back(d); // Start with the face of d Face f(v.dart); visitedFaces.push_back(f); // Start with the first face of v DartMarkerStore mark(map); mark.markOrbit(d) ; for(unsigned int iface = 0; iface != visitedFaces.size(); ++iface) mark.markOrbit(f) ; for(unsigned int iface = 0; iface != visitedFaces.size(); ++iface) { Dart e = visitedFaces[iface]; f = visitedFaces[iface]; VEC3 inter; bool interRes = intersectionLineConvexFace(map, e, position, point, dir, inter); bool interRes = intersectionLineConvexFace(map, f, position, point, dir, inter); if (interRes) { // check if already intersect on same point (a vertex certainly) bool alreadyfound = false; for(typename std::vector::iterator it = interPrec.begin(); !alreadyfound && it != interPrec.end(); ++it) { if (Geom::arePointsEquals(*it,inter)) alreadyfound = true; if (Geom::arePointsEquals(*it, inter)) alreadyfound = true; } if (!alreadyfound) { float v = dir * (inter-point); if (v>0) float v = dir * (inter - point); if (v > 0) ++countInter; if (v<0) if (v < 0) ++countInter2; interPrec.push_back(inter); } } // add all face neighbours to the table Dart currentFace = e; do foreach_adjacent2(map, f, [&] (Face ff) { Dart ee = map.phi2(e) ; if(!mark.isMarked(ee)) // not already marked if(!mark.isMarked(ff)) // not already marked { visitedFaces.push_back(ee) ; mark.markOrbit(ee) ; visitedFaces.push_back(ff) ; mark.markOrbit(ff) ; } e = map.phi1(e) ; } while(e != currentFace) ; }); } //if the point is in the volume there is an odd number of intersection with all faces with any direction return ((countInter % 2) != 0) && ((countInter2 % 2) != 0); // return (countInter % 2) == 1; } template bool isPointInConvexVolume(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const typename PFP::VEC3& point, bool CCW) bool isPointInConvexVolume(typename PFP::MAP& map, Vol v, const VertexAttribute& position, const typename PFP::VEC3& point, bool CCW) { typedef typename PFP::VEC3 VEC3 ; typedef typename PFP::REAL REAL; std::list visitedFaces; // Faces that are traversed visitedFaces.push_back(d); // Start with the face of d std::list::iterator face; VEC3 N; std::vector visitedFaces; // Faces that are traversed visitedFaces.reserve(64); Face f(v.dart); visitedFaces.push_back(f); // Start with the first face of v DartMarkerStore mark(map); // Lock a marker for (face = visitedFaces.begin(); face != visitedFaces.end(); ++face) for (std::vector::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face) { if (!mark.isMarked(*face)) f = *face; if (!mark.isMarked(f)) { Geom::Plane3D p = facePlane(map, *face, position); mark.markOrbit(f); Geom::Plane3D p = facePlane(map, f, position); Geom::Orientation3D o3d = p.orient(point); if(CCW) { ... ... @@ -152,15 +157,12 @@ bool isPointInConvexVolume(typename PFP::MAP& map, Dart d, const VertexAttribute else if(o3d == Geom::UNDER) return false; Dart dNext = *face ; do // add all face neighbours to the table foreach_adjacent2(map, f, [&] (Face ff) { mark.mark(dNext); // Mark Dart adj = map.phi2(dNext); // Get adjacent face if (adj != dNext && !mark.isMarked(adj)) visitedFaces.push_back(adj); // Add it dNext = map.phi1(dNext) ; } while(dNext != *face) ; if(!mark.isMarked(ff)) // not already marked visitedFaces.push_back(ff) ; }); } } ... ... @@ -168,22 +170,22 @@ bool isPointInConvexVolume(typename PFP::MAP& map, Dart d, const VertexAttribute } template bool isPointInConvexFace(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const typename PFP::VEC3& point, bool CCW) bool isPointInConvexFace(typename PFP::MAP& map, Face f, const VertexAttribute& position, const typename PFP::VEC3& point, bool CCW) { typedef typename PFP::VEC3 VEC3 ; typedef typename PFP::REAL REAL; Geom::Plane3D pl = Geometry::facePlane(map, d, position); Geom::Plane3D pl = Geometry::facePlane(map, f, position); Geom::Orientation3D o3d = pl.orient(point); if(o3d == Geom::ON) { Traversor2FV tfv(map, d) ; for(Dart it = tfv.begin(); it != tfv.end(); it = tfv.next()) Traversor2FV tfv(map, f) ; for(Vertex v = tfv.begin(); v != tfv.end(); v = tfv.next()) { VEC3 N = pl.normal(); VEC3 v2(position[map.phi1(it)] - position[it]); VEC3 norm2 = N^v2; Geom::Plane3D pl2(norm2, position[it]); VEC3 v2(position[map.phi1(v.dart)] - position[v]); VEC3 norm2 = N ^ v2; Geom::Plane3D pl2(norm2, position[v]); o3d = pl2.orient(point); if(CCW) { ... ... @@ -200,7 +202,7 @@ bool isPointInConvexFace(typename PFP::MAP& map, Dart d, const VertexAttribute bool isPointInConvexFace2D(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const typename PFP::VEC3& point, bool CCW ) bool isPointInConvexFace2D(typename PFP::MAP& map, Face f, const VertexAttribute& position, const typename PFP::VEC3& point, bool CCW ) { typedef typename PFP::VEC3 VEC3 ; typedef typename PFP::REAL REAL; ... ... @@ -209,10 +211,10 @@ bool isPointInConvexFace2D(typename PFP::MAP& map, Dart d, const VertexAttribute Geom::Orientation2D o2d; Traversor2FV tfv(map, d) ; for(Dart it = tfv.begin(); it != tfv.end(); it = tfv.next()) Traversor2FV tfv(map, f) ; for(Vertex v = tfv.begin(); v != tfv.end(); v = tfv.next()) { o2d = Geom::testOrientation2D(point, position[it], position[map.phi1(it)]); o2d = Geom::testOrientation2D(point, position[v], position[map.phi1(v.dart)]); if(CCW) { if(o2d == Geom::RIGHT) ... ... @@ -226,7 +228,7 @@ bool isPointInConvexFace2D(typename PFP::MAP& map, Dart d, const VertexAttribute } template bool isPointOnEdge(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const typename PFP::VEC3& point) bool isPointOnEdge(typename PFP::MAP& map, Edge e, const VertexAttribute& position, const typename PFP::VEC3& point) { // typedef typename PFP::REAL REAL; // typedef typename PFP::VEC3 VEC3 ; ... ... @@ -239,6 +241,8 @@ bool isPointOnEdge(typename PFP::MAP& map, Dart d, const VertexAttribute::min(); Dart d = e.dart; if( ( isPointOnHalfEdge(map,d,position,point) && isPointOnHalfEdge(map,map.phi2(d),position,point) ) || isPointOnVertex(map,d,position,point) || ... ... @@ -265,34 +269,34 @@ bool isPointOnHalfEdge(typename PFP::MAP& map, Dart d, const VertexAttribute= -REAL(0.00001); return abs(v1*v2) <= REAL(0.00001); } template bool isPointOnVertex(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const typename PFP::VEC3& point) bool isPointOnVertex(typename PFP::MAP& map, Vertex v, const VertexAttribute& position, const typename PFP::VEC3& point) { return Geom::arePointsEquals(point, position[d]); return Geom::arePointsEquals(point, position[v]); } template bool isConvexFaceInOrIntersectingTetrahedron(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const typename PFP::VEC3 points[4], bool CCW) bool isConvexFaceInOrIntersectingTetrahedron(typename PFP::MAP& map, Face f, const VertexAttribute& position, const typename PFP::VEC3 points[4], bool CCW) { typedef typename PFP::VEC3 VEC3 ; Traversor2FV tfv(map, d) ; for(Dart it = tfv.begin(); it != tfv.end(); it = tfv.next()) Traversor2FV tfv(map, f) ; for(Vertex v = tfv.begin(); v != tfv.end(); v = tfv.next()) { if(Geom::isPointInTetrahedron(points, position[it], CCW)) if(Geom::isPointInTetrahedron(points, position[v], CCW)) return true; } VEC3 inter; if( intersectionSegmentConvexFace(map, d, position, points[0], points[1], inter) || intersectionSegmentConvexFace(map, d, position, points[1], points[2], inter) || intersectionSegmentConvexFace(map, d, position, points[2], points[0], inter) || intersectionSegmentConvexFace(map, d, position, points[0], points[3], inter) || intersectionSegmentConvexFace(map, d, position, points[1], points[3], inter) || intersectionSegmentConvexFace(map, d, position, points[2], points[3], inter) if( intersectionSegmentConvexFace(map, f, position, points[0], points[1], inter) || intersectionSegmentConvexFace(map, f, position, points[1], points[2], inter) || intersectionSegmentConvexFace(map, f, position, points[2], points[0], inter) || intersectionSegmentConvexFace(map, f, position, points[0], points[3], inter) || intersectionSegmentConvexFace(map, f, position, points[1], points[3], inter) || intersectionSegmentConvexFace(map, f, position, points[2], points[3], inter) ) return true; ... ...
 ... ... @@ -42,47 +42,48 @@ namespace Geometry * (works with triangular faces but not optimized) * TODO to test * @param map the map * @param d a dart defining the n-sided face * @param f a n-sided face * @param PA segment point 1 * @param Dir a direction for the line * @param Inter store the intersection point * @return true if segment intersects the face */ template bool intersectionLineConvexFace(typename PFP::MAP& map, Dart d, const VertexAttribute& position, const typename PFP::VEC3& P, const typename PFP::VEC3& Dir, typename PFP::VEC3& Inter) ; bool intersectionLineConvexFace(typename PFP::MAP& map, Face f, const VertexAttribute& position, const typename PFP::VEC3& P, const typename PFP::VEC3& Dir, typename PFP::VEC3& Inter) ; /** * test the intersection between a segment and a n-sided face (n>=3) * (works with triangular faces but not optimized) * TODO optimize - based on intersectionLineConvexFace - check whether the points are on both sides of the face before computing intersection * @param map the map