First version of the cylinder clipping procedure.
Now it cut the mesh (the cap part is still lacking... See the app/trimesh_cylclip example
This commit is contained in:
parent
993a9a2c5b
commit
6762cdbacb
|
|
@ -0,0 +1,311 @@
|
|||
#ifndef CYLINDER_CLIP_H
|
||||
#define CYLINDER_CLIP_H
|
||||
#include <vcg/space/segment3.h>
|
||||
#include <vcg/complex/algorithms/refine.h>
|
||||
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-<v,p1-p+td1>d)^2 -r^2=0, becomes At^2+Bt+C=0.
|
||||
//
|
||||
// tmpA = d1 - (<d1,d>/<d,d>)*d.
|
||||
// tmpB = (p1-p) - (<p1-p,d>/<d,d>)*d.
|
||||
// A = <tmpA,tmpA>.
|
||||
// B = 2*<tmpA,tmpB>.
|
||||
// C = <tmpB,tmpB> - r^2.
|
||||
// Input: Cylinder<T> & cyl, Line3<T> & line.
|
||||
// Output: CoordType & p0,CoordType & p1, T & t0, T &t1.
|
||||
|
||||
template<class T>
|
||||
static bool IntersectionLineCylinder(const Segment3<T> & cylSeg, T radius, const Line3<T> & line, Point3<T> & p0, Point3<T> & p1, T & t0, T &t1)
|
||||
{
|
||||
T dist;
|
||||
Point3<T> mClosestPoint0,mClosestPoint1;
|
||||
bool parallel=false;
|
||||
|
||||
Line3<T> 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<T> cyl_origin=tmp.Origin();
|
||||
Point3<T> line_origin=line.Origin();
|
||||
Point3<T> cyl_direction=tmp.Direction();
|
||||
Point3<T> line_direction=line.Direction();
|
||||
|
||||
Point3<T> deltaP=line_origin-cyl_origin;
|
||||
|
||||
T dotDirCyl=cyl_direction.SquaredNorm(); //<d,d>
|
||||
|
||||
T scalar=line_direction.dot(cyl_direction);
|
||||
Point3<T> tmpA=line_direction-(cyl_direction/dotDirCyl)*scalar;
|
||||
T A=tmpA.SquaredNorm();
|
||||
|
||||
T scalar2=deltaP.dot(cyl_direction);
|
||||
Point3<T> 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<MESH_TYPE, Type, T> & cyl, Segment3<T> & seg.
|
||||
// Output: CoordType & p0,CoordType & p1, T & t0, T &t1.
|
||||
|
||||
template<class T>
|
||||
static bool IntersectionSegmentCylinder(const Segment3<T> & cylSeg , T radius, const Segment3<T> & seg, Point3<T> & p0, Point3<T> & p1)
|
||||
{
|
||||
const float eps = 10e-5;
|
||||
|
||||
Line3<T> 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<class T>
|
||||
static bool PointIsInSegment(const Point3<T> &point, const Segment3<T> &seg){
|
||||
const float eps = 10e-5;
|
||||
Line3<T> 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 MeshType>
|
||||
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<FaceType> PosType;
|
||||
typedef Segment3<ScalarType> Segment3x;
|
||||
typedef Plane3<ScalarType> 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,CoordType>,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??,
|
||||
// is external if the distance is > radius+eps and it is on the circumference if the
|
||||
// distance is in the range [radius-eps??, radius+eps].
|
||||
//
|
||||
// Input: face::Pos<typename MESH_TYPE::FaceType> ep, Cylinder<typename MESH_TYPE::ScalarType> 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)<maxDist && fabs(dist1-radius)<maxDist)
|
||||
{
|
||||
newPtMap[std::make_pair(p0,p1)] = (p0+p1)*0.5;
|
||||
return true;
|
||||
}
|
||||
|
||||
// ************ Case 1;
|
||||
if((dist0>radius) && (dist1>radius))
|
||||
{
|
||||
bool parallel;
|
||||
SegmentSegmentDistance(edgeSeg,this->seg, dist2, parallel, closest0,closest1);
|
||||
if((dist2<radius) &&
|
||||
(Distance(closest0,p0)>minEdgeLen) &&
|
||||
(Distance(closest0,p1)>minEdgeLen))
|
||||
{
|
||||
newPtMap[std::make_pair(p0,p1)] = closest0;
|
||||
return true;
|
||||
}
|
||||
}
|
||||
else if(((dist0<radius) && (dist1>radius))||((dist0>radius) && (dist1<radius))){
|
||||
CoordType int0,int1;
|
||||
// If there is an intersection point between segment and cylinder,
|
||||
// this must be different from the extremes of the segment and
|
||||
// his projection must be in the segment.
|
||||
if(IntersectionSegmentCylinder(this->seg, this->radius,edgeSeg,int0,int1)){
|
||||
if(PointIsInSegment(int0,this->seg) && (Distance(p0,int0)>eps) && (Distance(p1,int0)>eps))
|
||||
{
|
||||
if(Distance(int0,p0)<maxDist) return false;
|
||||
if(Distance(int0,p1)<maxDist) return false;
|
||||
newPtMap[std::make_pair(p0,p1)] = int0;
|
||||
return true;
|
||||
}
|
||||
}
|
||||
}
|
||||
// Now check also against the caps
|
||||
CoordType pt;
|
||||
if(IntersectionPlaneSegment(pl0,edgeSeg,pt)){
|
||||
if((Distance(pt,origin)<radius+1.0*minEdgeLen) &&
|
||||
(Distance(pt,p0)>eps) && (Distance(pt,p1)>eps) )
|
||||
{
|
||||
newPtMap[std::make_pair(p0,p1)] = pt;
|
||||
return true;
|
||||
}
|
||||
}
|
||||
if(IntersectionPlaneSegment(pl1,edgeSeg,pt)){
|
||||
if( (Distance(pt,end)<radius+1.0*minEdgeLen) &&
|
||||
(Distance(pt,p0)>eps) && (Distance(pt,p1)>eps) )
|
||||
{
|
||||
newPtMap[std::make_pair(p0,p1)] = pt;
|
||||
return true;
|
||||
}
|
||||
}
|
||||
return false;
|
||||
//
|
||||
}
|
||||
};
|
||||
|
||||
class CylMidPoint : public std::unary_function<PosType, CoordType>
|
||||
{
|
||||
private:
|
||||
CylMidPoint() {assert(0);}
|
||||
public:
|
||||
CylMidPoint(CylPred &ep) : newPtMap(&(ep.newPtMap)) {
|
||||
assert(newPtMap);
|
||||
}
|
||||
std::map< std::pair<CoordType,CoordType>, CoordType > *newPtMap;
|
||||
void operator()(VertexType &nv, PosType ep)
|
||||
{
|
||||
typename std::map< std::pair<CoordType,CoordType>,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<ScalarType> WedgeInterp(Color4<ScalarType> &c0, Color4<ScalarType> &c1)
|
||||
{
|
||||
Color4<ScalarType> cc;
|
||||
return cc.lerp(c0,c1,0.5f);
|
||||
}
|
||||
|
||||
TexCoord2<ScalarType,1> WedgeInterp(TexCoord2<ScalarType,1> &t0, TexCoord2<ScalarType,1> &t1)
|
||||
{
|
||||
TexCoord2<ScalarType,1> 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<MeshType, CylMidPoint >(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
|
||||
Loading…
Reference in New Issue