From 6762cdbacb8aec5f1c74002f2250a07005977528 Mon Sep 17 00:00:00 2001 From: cignoni Date: Thu, 28 Nov 2013 23:30:35 +0000 Subject: [PATCH] First version of the cylinder clipping procedure. Now it cut the mesh (the cap part is still lacking... See the app/trimesh_cylclip example --- vcg/complex/algorithms/cylinder_clipping.h | 311 +++++++++++++++++++++ 1 file changed, 311 insertions(+) create mode 100644 vcg/complex/algorithms/cylinder_clipping.h diff --git a/vcg/complex/algorithms/cylinder_clipping.h b/vcg/complex/algorithms/cylinder_clipping.h new file mode 100644 index 00000000..241ea1c5 --- /dev/null +++ b/vcg/complex/algorithms/cylinder_clipping.h @@ -0,0 +1,311 @@ +#ifndef CYLINDER_CLIP_H +#define CYLINDER_CLIP_H +#include +#include +namespace vcg +{ + +// Taken a cylinder and a line calculates whether there is an intersection between these. +// In output are provided, if they exist, any two points of intersection (p0, p1) +// and the parameters t (t0, t1) on the line. +// Returns false if the distance of the line from the axis of the cylinder is greater than +// the radius of the cylinder or, if the calculation of t parameters - obtained by solving the +// quadratic equation - gives a delta less than zero. +// To find the intersection of a line p1+td1 with the axis p+td of the cylinder: +// (p1-p+td1-d)^2 -r^2=0, becomes At^2+Bt+C=0. +// +// tmpA = d1 - (/)*d. +// tmpB = (p1-p) - (/)*d. +// A = . +// B = 2*. +// C = - r^2. +// Input: Cylinder & cyl, Line3 & line. +// Output: CoordType & p0,CoordType & p1, T & t0, T &t1. + +template +static bool IntersectionLineCylinder(const Segment3 & cylSeg, T radius, const Line3 & line, Point3 & p0, Point3 & p1, T & t0, T &t1) +{ + T dist; + Point3 mClosestPoint0,mClosestPoint1; + bool parallel=false; + + Line3 tmp; + tmp.Set(cylSeg.P0(), (cylSeg.P1()-cylSeg.P0()).Normalize()); + LineLineDistance(tmp,line,parallel,dist,mClosestPoint0,mClosestPoint1); + if(dist>radius) + return false; + if(parallel) return false; + Point3 cyl_origin=tmp.Origin(); + Point3 line_origin=line.Origin(); + Point3 cyl_direction=tmp.Direction(); + Point3 line_direction=line.Direction(); + + Point3 deltaP=line_origin-cyl_origin; + + T dotDirCyl=cyl_direction.SquaredNorm(); // + + T scalar=line_direction.dot(cyl_direction); + Point3 tmpA=line_direction-(cyl_direction/dotDirCyl)*scalar; + T A=tmpA.SquaredNorm(); + + T scalar2=deltaP.dot(cyl_direction); + Point3 tmpB=(deltaP-(cyl_direction/dotDirCyl)*scalar2); + + T B=2.0*tmpA.dot(tmpB); + + T C=tmpB.SquaredNorm()-pow(radius,2); + + T delta=pow(B,2)-4*A*C; + + if(delta<0) + return false; + + t0=(-B-sqrt(delta))/(2*A); + t1=(-B+sqrt(delta))/(2*A); + + p0=line.P(t0); + p1=line.P(t1); + + return true; +} + +// Taken a cylinder and a segment calculates the intersection possible using the +// IntersectionLineCylinder() and checking the output of this. +// Whether the t0 and t1 scalars are between 0 and the length of the segment, then the point +// belongs to it and returns true. +// In output are given two points of intersection (p0, p1) and the parameters t (t0, t1) on the line. +// If p1 belongs to the segment and p0 no, it swaps the points (p0, p1) because operator() in the +// MidPointCylinder always takes the first. +// Otherwise, it means that there is no point between the extremes of the segment that intersects +// the cylinder, in this case it returns false. +// +// Input: Cylinder & cyl, Segment3 & seg. +// Output: CoordType & p0,CoordType & p1, T & t0, T &t1. + +template +static bool IntersectionSegmentCylinder(const Segment3 & cylSeg , T radius, const Segment3 & seg, Point3 & p0, Point3 & p1) +{ + const float eps = 10e-5; + + Line3 line; + line.SetOrigin(seg.P0()); + line.SetDirection((seg.P1()-seg.P0()).Normalize()); + + T t0,t1; + if(IntersectionLineCylinder(cylSeg,radius,line,p0,p1,t0,t1)){ + bool inters0 = (t0>=0) && (t0<=seg.Length()); + bool inters1 = (t1>=0) && (t1<=seg.Length()); + if( inters0 && !inters1) p1=p0; // if only one of the line intersections belong to the segment + if(!inters0 && inters1) p0=p1; // make both of them the same value. + return inters0 || inters1; + } + return false; +} + +template +static bool PointIsInSegment(const Point3 &point, const Segment3 &seg){ + const float eps = 10e-5; + Line3 line; + line.SetOrigin(seg.P0()); + line.SetDirection((seg.P1()-seg.P0())); + T t=line.Projection(point); + // Remembers, two points are different if their distance is >=eps + if(t>-eps && t<1+eps) + return true; + return false; +} + +namespace tri +{ + +template +class CylinderClipping +{ +public: + typedef typename MeshType::ScalarType ScalarType; + typedef typename MeshType::VertexType VertexType; + typedef typename MeshType::VertexPointer VertexPointer; + typedef typename MeshType::VertexIterator VertexIterator; + typedef typename MeshType::CoordType CoordType; + typedef typename MeshType::FaceType FaceType; + typedef typename MeshType::FacePointer FacePointer; + typedef typename MeshType::FaceIterator FaceIterator; + typedef typename face::Pos PosType; + typedef Segment3 Segment3x; + typedef Plane3 Plane3x; + + // This predicate + class CylPred + { + public: + CylPred(CoordType &_origin, CoordType &_end, ScalarType _radius, ScalarType _maxDist, ScalarType _minEdgeLen): + origin(_origin),end(_end),radius(_radius),maxDist(_maxDist),minEdgeLen(_minEdgeLen){ + seg.Set(origin,end); + pl0.Init(origin,(end-origin).Normalize()); + pl1.Init(end,(end-origin).Normalize()); + } + void Init() { newPtMap.clear(); } + ScalarType radius; + CoordType origin,end; + ScalarType minEdgeLen; + ScalarType maxDist; + private: + Segment3x seg; + Plane3x pl0,pl1; + public: + // This map store for each edge the new position. + // it is intializaed by the predicate itself. + // and it is shared with the midpoint functor. + std::map< std::pair,CoordType > newPtMap; + + // Return true if the given edge intersect the cylinder. + + // Verify if exist a point in an edge that intersects the cylinder. Then calculate + // this point and store it for later use. + // The cases possible are: + // 1. Both extremes have distance greater than or equal to the radius, in this case it + // calculates the point of this segment closest to the axis of the cylinder. If this + // has distance less than or equal to the radius and is different from the extremes + // returns true and this point, otherwise false; + // 2. If there is an extreme inside and one outside it returns true because exist the point + // of intersection that is calculated using the IntersectionSegmentCylinder(); + // 3. Otherwise false. + // So a point is inside of the cylinder if its distance from his axis is radius+eps and it is on the circumference if the + // distance is in the range [radius-eps??, radius+eps]. + // + // Input: face::Pos ep, Cylinder cyl, + + bool operator()(PosType ep) + { + VertexType *v0 = ep.V(); + VertexType *v1 = ep.VFlip(); + ScalarType eps = minEdgeLen/100.0f; + + if(v0>v1) std::swap(v0,v1); + CoordType p0=v0->P(); + CoordType p1=v1->P(); + + // CASE 0 - For very short edges --- DO NOTHING + if(Distance(p0,p1)< minEdgeLen) return false; + + Segment3x edgeSeg(p0,p1); + CoordType closest0,closest1; // points on the cyl axis + ScalarType dist0,dist1,dist2; + + SegmentPointDistance(this->seg,p0,closest0,dist0); + SegmentPointDistance(this->seg,p1,closest1,dist1); + + // Case 0.5 + if(fabs(dist0-radius)radius) && (dist1>radius)) + { + bool parallel; + SegmentSegmentDistance(edgeSeg,this->seg, dist2, parallel, closest0,closest1); + if((dist2minEdgeLen) && + (Distance(closest0,p1)>minEdgeLen)) + { + newPtMap[std::make_pair(p0,p1)] = closest0; + return true; + } + } + else if(((dist0radius))||((dist0>radius) && (dist1seg, this->radius,edgeSeg,int0,int1)){ + if(PointIsInSegment(int0,this->seg) && (Distance(p0,int0)>eps) && (Distance(p1,int0)>eps)) + { + if(Distance(int0,p0)eps) && (Distance(pt,p1)>eps) ) + { + newPtMap[std::make_pair(p0,p1)] = pt; + return true; + } + } + if(IntersectionPlaneSegment(pl1,edgeSeg,pt)){ + if( (Distance(pt,end)eps) && (Distance(pt,p1)>eps) ) + { + newPtMap[std::make_pair(p0,p1)] = pt; + return true; + } + } + return false; + // + } + }; + + class CylMidPoint : public std::unary_function + { + private: + CylMidPoint() {assert(0);} + public: + CylMidPoint(CylPred &ep) : newPtMap(&(ep.newPtMap)) { + assert(newPtMap); + } + std::map< std::pair, CoordType > *newPtMap; + void operator()(VertexType &nv, PosType ep) + { + typename std::map< std::pair,CoordType >::iterator mi; + VertexType *v0 = ep.V(); + VertexType *v1 = ep.VFlip(); + assert(newPtMap); + if(v0>v1) std::swap(v0,v1); + + CoordType p0=v0->P(); + CoordType p1=v1->P(); + mi=newPtMap->find(std::make_pair(v0->P(),v1->P())); + assert(mi!=newPtMap->end()); + nv.P()=(*mi).second; + } + + Color4 WedgeInterp(Color4 &c0, Color4 &c1) + { + Color4 cc; + return cc.lerp(c0,c1,0.5f); + } + + TexCoord2 WedgeInterp(TexCoord2 &t0, TexCoord2 &t1) + { + TexCoord2 tmp; + assert(t0.n()== t1.n()); + tmp.n()=t0.n(); + tmp.t()=(t0.t()+t1.t())/2.0; + return tmp; + } + }; + + static void Apply(MeshType &m, CoordType &origin, CoordType &end, ScalarType radius) + { + CylPred cylep(origin,end,radius,radius/100.0,m.bbox.Diag()/50.0f); + CylMidPoint cylmp(cylep); + int i=0; + while((tri::RefineE(m, cylmp,cylep))&&(i<50)){ + cylep.Init(); + printf("Refine %d Vertici: %d, Facce: %d\n",i,m.VN(),m.FN()); + i++; + } + } +}; +} // end namespace tri +} // end namespace vcg +#endif // CYLINDER_CLIP_H