Commit fee3b1fc by Pierre Kraemer

### corrections intersection sphere / edge

parent b18a529d
 ... ... @@ -360,8 +360,9 @@ void computeCurvatureVertex_NormalCycles( const std::vector& vd2 = neigh.getBorder() ; for (std::vector::const_iterator it = vd2.begin(); it != vd2.end(); ++it) { const VEC3 e = position[map.phi2(*it)] - position[*it] ; const REAL alpha = neigh.intersect_SphereEdge(*it, map.phi2(*it)) ; const VEC3 e = Algo::Geometry::vectorOutOfDart(map, *it, position) ; REAL alpha ; Algo::Geometry::intersectionSphereEdge(map, center, radius, *it, position, alpha) ; tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) * alpha ; } ... ...
 ... ... @@ -50,12 +50,12 @@ namespace Geometry * @return true if segment intersects the face */ template bool intersectionLineConvexFace(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& positions, const typename PFP::VEC3& PA, const typename PFP::VEC3& Dir, typename PFP::VEC3& Inter) ; bool intersectionLineConvexFace(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& 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 to optimize -based on intersectionLineConvexFace- check whether the points are on both sides of the face before computing intersection * TODO optimize - based on intersectionLineConvexFace - check whether the points are on both sides of the face before computing intersection * @param map the map * @param d a dart defining the n-sided face * @param PA segment point 1 ... ... @@ -64,7 +64,7 @@ bool intersectionLineConvexFace(typename PFP::MAP& map, Dart d, const typename P * @return true if segment intersects the face */ template bool intersectionSegmentConvexFace(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& positions, const typename PFP::VEC3& PA, const typename PFP::VEC3& PB, typename PFP::VEC3& Inter) ; bool intersectionSegmentConvexFace(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position, const typename PFP::VEC3& PA, const typename PFP::VEC3& PB, typename PFP::VEC3& Inter) ; /** * test if two triangles intersect ... ... @@ -74,13 +74,16 @@ bool intersectionSegmentConvexFace(typename PFP::MAP& map, Dart d, const typenam * @param a dart of the second triangle */ template bool areTrianglesInIntersection(typename PFP::MAP& map, Dart tri1, Dart tri2, const typename PFP::TVEC3& positions) ; bool areTrianglesInIntersection(typename PFP::MAP& map, Dart tri1, Dart tri2, const typename PFP::TVEC3& position) ; /** * * alpha = coef d'interpolation dans [0 ,1] tel que v = (1-alpha)*pin + alpha*pout * est le point d'intersection entre la sphère et le segment [pin, pout] * avec pin = position[din] à l'intérieur de la sphère * avec pout = position[dout] à l'extérieur de la sphère */ template bool intersectionSphereEdge(typename PFP::MAP& map, typename PFP::VEC3& center, typename PFP::REAL radius, Dart d, const typename PFP::TVEC3& positions, typename PFP::REAL& alpha) ; bool intersectionSphereEdge(typename PFP::MAP& map, typename PFP::VEC3& center, typename PFP::REAL radius, Dart d, const typename PFP::TVEC3& position, typename PFP::REAL& alpha) ; } // namespace Geometry ... ...
 ... ... @@ -24,7 +24,6 @@ #include "Algo/Geometry/normal.h" #include "Algo/Geometry/centroid.h" #include "intersection.h" #include ... ... @@ -38,15 +37,14 @@ namespace Geometry { template bool intersectionLineConvexFace(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& positions, const typename PFP::VEC3& P, const typename PFP::VEC3& Dir, typename PFP::VEC3& Inter) bool intersectionLineConvexFace(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position, const typename PFP::VEC3& P, const typename PFP::VEC3& Dir, typename PFP::VEC3& Inter) { typedef typename PFP::VEC3 VEC3 ; const float SMALL_NUM = std::numeric_limits::min() * 5.0f; VEC3 p1 = positions[d]; VEC3 n = faceNormal(map, d, positions); VEC3 p1 = position[d]; VEC3 n = faceNormal(map, d, position); VEC3 w0 = P - p1; float a = -(n * w0); float b = n * Dir; ... ... @@ -58,7 +56,7 @@ bool intersectionLineConvexFace(typename PFP::MAP& map, Dart d, const typename P Inter = P + r * Dir; // intersect point of ray and plane // is I inside the face? VEC3 p2 = positions[map.phi1(d)]; VEC3 p2 = position[map.phi1(d)]; VEC3 v = p2 - p1 ; VEC3 vInter = Inter - p1; float dirV = v * vInter; ... ... @@ -69,7 +67,7 @@ bool intersectionLineConvexFace(typename PFP::MAP& map, Dart d, const typename P while(it != d) { p1 = p2; p2 = positions[map.phi1(it)]; p2 = position[map.phi1(it)]; v = p2 - p1; vInter = Inter - p1; float dirD = v * vInter; ... ... @@ -85,22 +83,22 @@ bool intersectionLineConvexFace(typename PFP::MAP& map, Dart d, const typename P } template bool intersectionSegmentConvexFace(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& positions, const typename PFP::VEC3& PA, const typename PFP::VEC3& PB, typename PFP::VEC3& Inter) bool intersectionSegmentConvexFace(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position, const typename PFP::VEC3& PA, const typename PFP::VEC3& PB, typename PFP::VEC3& Inter) { typename PFP::VEC3 dir = PB - PA; if (intersectionLineConvexFace(map,d,positions,PA,dir,Inter)) if (intersectionLineConvexFace(map, d, position, PA, dir, Inter)) { typename PFP::VEC3 dirA = PA - Inter; typename PFP::VEC3 dirB = PB -Inter; typename PFP::VEC3 dirB = PB - Inter; if ( (dirA*dirB) < 0 ) if (dirA * dirB < 0) return true; } return false; } template bool areTrianglesInIntersection(typename PFP::MAP& map, Dart tri1, Dart tri2, const typename PFP::TVEC3& positions) bool areTrianglesInIntersection(typename PFP::MAP& map, Dart tri1, Dart tri2, const typename PFP::TVEC3& position) { typedef typename PFP::VEC3 VEC3 ; ... ... @@ -109,8 +107,8 @@ bool areTrianglesInIntersection(typename PFP::MAP& map, Dart tri1, Dart tri2, co VEC3 tris2[3]; for (unsigned int i = 0; i < 3; ++i) { tris1[i] = positions[tri1]; tris2[i] = positions[tri2]; tris1[i] = position[tri1]; tris2[i] = position[tri2]; tri1 = map.phi1(tri1); tri2 = map.phi1(tri2); } ... ... @@ -167,8 +165,8 @@ bool areTrianglesInIntersection(typename PFP::MAP& map, Dart tri1, Dart tri2, co // return isSegmentInTriangle2D(inter1,inter2,tris2[0],tris2[1],triS2[2],nTri2); //compute face normal VEC3 normale1 = faceNormal(map,tri1,positions); VEC3 bary1 = faceCentroid(map,tri1,positions); VEC3 normale1 = faceNormal(map, tri1, position); VEC3 bary1 = faceCentroid(map, tri1, position); int pos = 0; int neg = 0; ... ... @@ -188,8 +186,8 @@ bool areTrianglesInIntersection(typename PFP::MAP& map, Dart tri1, Dart tri2, co return false; //same for the second triangle VEC3 normale2 = faceNormal(map,tri2,positions); VEC3 bary2 = faceCentroid(map,tri2,positions); VEC3 normale2 = faceNormal(map, tri2, position); VEC3 bary2 = faceCentroid(map, tri2, position); pos = 0; neg = 0; for (unsigned int i = 0; i < 3 ; ++i) ... ... @@ -210,7 +208,7 @@ bool areTrianglesInIntersection(typename PFP::MAP& map, Dart tri1, Dart tri2, co for (unsigned int i = 0; i < 3 && !intersection; ++i) { VEC3 inter; intersection = Geom::intersectionSegmentTriangle(tris1[i],tris1[(i+1)%3],tris2[0],tris2[1],tris2[2],inter); intersection = Geom::intersectionSegmentTriangle(tris1[i], tris1[(i+1)%3], tris2[0], tris2[1], tris2[2], inter); } if (intersection) ... ... @@ -219,17 +217,20 @@ bool areTrianglesInIntersection(typename PFP::MAP& map, Dart tri1, Dart tri2, co for (unsigned int i = 0; i < 3 && !intersection; ++i) { VEC3 inter; intersection = Geom::intersectionSegmentTriangle(tris2[i],tris2[(i+1)%3],tris1[0],tris1[1],tris1[2],inter); intersection = Geom::intersectionSegmentTriangle(tris2[i], tris2[(i+1)%3], tris1[0], tris1[1], tris1[2], inter); } return intersection; } template bool intersectionSphereEdge(typename PFP::MAP& map, typename PFP::VEC3& center, typename PFP::REAL radius, Dart d, const typename PFP::TVEC3& positions, typename PFP::REAL& alpha) bool intersectionSphereEdge(typename PFP::MAP& map, typename PFP::VEC3& center, typename PFP::REAL radius, Dart d, const typename PFP::TVEC3& position, typename PFP::REAL& alpha) { typename PFP::VEC3& p1 = position[d]; typename PFP::VEC3& p2 = position[map.phi1(d)]; typedef typename PFP::VEC3 VEC3 ; typedef typename PFP::REAL REAL ; const typename PFP::VEC3& p1 = position[d]; const typename PFP::VEC3& p2 = position[map.phi1(d)]; if(Geom::isPointInSphere(p1, center, radius) && !Geom::isPointInSphere(p2, center, radius)) { VEC3 p = p1 - center; ... ...
 ... ... @@ -25,6 +25,8 @@ #ifndef __COLLECTOR_H__ #define __COLLECTOR_H__ #include "Algo/Geometry/intersection.h" /***************************************** * Class hierarchy : * Collector (virtual) ... ... @@ -135,7 +137,6 @@ template class Collector_WithinSphere : public Collector { protected: typename PFP::REAL radius_2; typename PFP::VEC3 centerPosition; typename PFP::REAL area; ... ... @@ -145,16 +146,6 @@ public: {} void init(Dart d, typename PFP::REAL r = 0); void collect(); bool isInside(Dart d) { return (this->position[d] - centerPosition).norm2() < radius_2; } // alpha = coef d'interpolation dans [0 ,1] tel que v= (1-alpha)*pin + alpha*pout // est le point d'intersection entre la sphère et le segment [pin, pout] // avec pin = position[din] à l'intérieur de la sphère // avec pout = position[dout] à l'extérieur de la sphère typename PFP::REAL intersect_SphereEdge(const Dart din, const Dart dout); void computeArea(); typename PFP::REAL getArea() const { return area; } }; ... ...
 ... ... @@ -112,7 +112,6 @@ void Collector_WithinSphere::init(Dart d, typename PFP::REAL r) this->insideEdges.clear(); this->insideFaces.clear(); this->border.clear(); radius_2 = r * r; centerPosition = this->position[d]; area = 0; } ... ... @@ -143,7 +142,7 @@ void Collector_WithinSphere::collect() const Dart f = this->map.phi1(e); const Dart g = this->map.phi1(f); if (! this->isInside(f)) if (! Geom::isPointInSphere(this->position[f], centerPosition, this->radius)) { this->border.push_back(e); // add to border em.mark(e); ... ... @@ -161,7 +160,7 @@ void Collector_WithinSphere::collect() this->insideEdges.push_back(e); em.mark(e); } if (! fm.isMarked(e) && this->isInside(g)) if (! fm.isMarked(e) && Geom::isPointInSphere(this->position[g], centerPosition, this->radius)) { this->insideFaces.push_back(e); fm.mark(e); ... ... @@ -174,41 +173,31 @@ void Collector_WithinSphere::collect() } } //template //typename PFP::REAL Collector_WithinSphere::intersect_SphereEdge(const Dart din, const Dart dout) //{ // typedef typename PFP::VEC3 VEC3; // typedef typename PFP::REAL REAL; // VEC3 p = this->position[din] - this->centerPosition; // VEC3 qminusp = this->position[dout] - this->centerPosition - p; // REAL s = p * qminusp; // REAL n2 = qminusp.norm2(); // return (- s + sqrt(s * s + n2 * (this->radius_2 - p.norm2()))) / n2; //} template void Collector_WithinSphere::computeArea() { for (std::vector::const_iterator it = this->insideFaces.begin(); it != this->insideFaces.end(); ++it) { this->area += Algo::Geometry::triangleArea(this->map, *it, this->position); area += Algo::Geometry::triangleArea(this->map, *it, this->position); } for (std::vector::const_iterator it = this->border.begin(); it != this->border.end(); ++it) { const Dart f = this->map.phi1(*it); // we know that f is outside const Dart g = this->map.phi1(f); if (this->isInside(g)) if (Geom::isPointInSphere(this->position[g], centerPosition, this->radius)) { // only f is outside typename PFP::REAL alpha = this->intersect_SphereEdge(*it, f); typename PFP::REAL beta = this->intersect_SphereEdge(g, f); this->area += (alpha+beta - alpha*beta) * Algo::Geometry::triangleArea(this->map, *it, this->position); typename PFP::REAL alpha, beta; Algo::Geometry::intersectionSphereEdge(this->map, centerPosition, this->radius, *it, this->position, alpha); Algo::Geometry::intersectionSphereEdge(this->map, centerPosition, this->radius, this->map.phi2(f), this->position, beta); area += (alpha+beta - alpha*beta) * Algo::Geometry::triangleArea(this->map, *it, this->position); } else { // f and g are outside typename PFP::REAL alpha = this->intersect_SphereEdge(*it, f); typename PFP::REAL beta = this->intersect_SphereEdge(*it, g); this->area += alpha * beta * Algo::Geometry::triangleArea(this->map, *it, this->position); typename PFP::REAL alpha, beta; Algo::Geometry::intersectionSphereEdge(this->map, centerPosition, this->radius, *it, this->position, alpha); Algo::Geometry::intersectionSphereEdge(this->map, centerPosition, this->radius, this->map.phi2(g), this->position, beta); area += alpha * beta * Algo::Geometry::triangleArea(this->map, *it, this->position); } } } ... ...
 ... ... @@ -54,7 +54,7 @@ template Inclusion isPointInTriangle(const VEC3& point, const VEC3& Ta, const VEC3& Tb, const VEC3& Tc) ; template Inclusion isPointInSphere(const VEC3& point, const VEC3& center, const typename VEC3::DATA_TYPE& radius) ; bool isPointInSphere(const VEC3& point, const VEC3& center, const typename VEC3::DATA_TYPE& radius) ; /** * test if a segment is inside a triangle, the segment MUST be in the plane of the triangle ... ...
 ... ... @@ -95,9 +95,9 @@ Intersection intersectionPlaneRay(const PLANE3D& pl,const VEC3& p1,const VEC3& d template Intersection intersection2DSegmentSegment(const VEC3& PA, const VEC3& PB, const VEC3& QA, const VEC3& QB, VEC3& Inter) ; } } // namespace Geom } } // namespace CGoGN #include "intersection.hpp" ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!