+ /*!
+ * \brief Compute point position on a curve. Use octree to fast reject far points
+ */
+ class CurveProjector : public SMESH_Octree
+ {
+ public:
+ CurveProjector( const TopoDS_Edge& edge, double enlarge );
+
+ bool IsOnCurve( const gp_XYZ& point, double & distance2, double & u );
+
+ bool IsOut( const gp_XYZ& point ) const { return getBox()->IsOut( point ); }
+
+ protected:
+ CurveProjector() {}
+ SMESH_Octree* newChild() const { return new CurveProjector; }
+ void buildChildrenData();
+ Bnd_B3d* buildRootBox();
+
+ private:
+ struct CurveSegment : public Bnd_B3d
+ {
+ double _chord, _chord2, _length2;
+ gp_Pnt _pFirst, _pLast;
+ gp_Lin _line;
+ Handle(Geom_Curve) _curve;
+
+ CurveSegment() {}
+ void Init( const gp_Pnt& pf, const gp_Pnt& pl,
+ double uf, double ul, double tol, Handle(Geom_Curve)& curve );
+ bool IsOn( const gp_XYZ& point, double & distance2, double & u );
+ bool IsInContact( const Bnd_B3d& bb );
+ };
+ std::vector< CurveSegment > _segments;
+ };
+
+ //===============================================================================
+ /*!
+ * \brief Create an octree of curve segments
+ */
+ //================================================================================
+
+ CurveProjector::CurveProjector( const TopoDS_Edge& edge, double enlarge )
+ :SMESH_Octree( 0 )
+ {
+ double f,l;
+ Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, f, l );
+ double curDeflect = 0.3; // Curvature deflection
+ double angDeflect = 1e+100; // Angular deflection - don't control chordal error
+ GCPnts_TangentialDeflection div( BRepAdaptor_Curve( edge ), angDeflect, curDeflect );
+ _segments.resize( div.NbPoints() - 1 );
+ for ( int i = 1; i < div.NbPoints(); ++i )
+ try {
+ _segments[ i - 1 ].Init( div.Value( i ), div.Value( i+1 ),
+ div.Parameter( i ), div.Parameter( i+1 ),
+ enlarge, curve );
+ }
+ catch ( Standard_Failure ) {
+ _segments.resize( _segments.size() - 1 );
+ --i;
+ }
+ if ( _segments.size() < 3 )
+ myIsLeaf = true;
+
+ compute();
+
+ if ( _segments.size() == 1 )
+ myBox->Enlarge( enlarge );
+ }
+
+ //================================================================================
+ /*!
+ * \brief Return the maximal box
+ */
+ //================================================================================
+
+ Bnd_B3d* CurveProjector::buildRootBox()
+ {
+ Bnd_B3d* box = new Bnd_B3d;
+ for ( size_t i = 0; i < _segments.size(); ++i )
+ box->Add( _segments[i] );
+ return box;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Redistribute segments among children
+ */
+ //================================================================================
+
+ void CurveProjector::buildChildrenData()
+ {
+ bool allIn = true;
+ for ( size_t i = 0; i < _segments.size(); ++i )
+ {
+ for (int j = 0; j < 8; j++)
+ {
+ if ( _segments[i].IsInContact( *myChildren[j]->getBox() ))
+ ((CurveProjector*)myChildren[j])->_segments.push_back( _segments[i]);
+ else
+ allIn = false;
+ }
+ }
+ if ( allIn && _segments.size() < 3 )
+ {
+ myIsLeaf = true;
+ for (int j = 0; j < 8; j++)
+ static_cast<CurveProjector*>( myChildren[j])->myIsLeaf = true;
+ }
+ else
+ {
+ SMESHUtils::FreeVector( _segments ); // = _segments.clear() + free memory
+
+ for (int j = 0; j < 8; j++)
+ {
+ CurveProjector* child = static_cast<CurveProjector*>( myChildren[j]);
+ if ( child->_segments.size() < 3 )
+ child->myIsLeaf = true;
+ }
+ }
+ }
+
+ //================================================================================
+ /*!
+ * \brief Return true if a point is close to the curve
+ * \param [in] point - the point
+ * \param [out] distance2 - distance to the curve
+ * \param [out] u - parameter on the curve
+ * \return bool - is the point is close to the curve
+ */
+ //================================================================================
+
+ bool CurveProjector::IsOnCurve( const gp_XYZ& point, double & distance2, double & u )
+ {
+ if ( getBox()->IsOut( point ))
+ return false;
+
+ bool ok = false;
+ double dist2, param;
+ distance2 = Precision::Infinite();
+
+ if ( isLeaf() )
+ {
+ for ( size_t i = 0; i < _segments.size(); ++i )
+ if ( !_segments[i].IsOut( point ) &&
+ _segments[i].IsOn( point, dist2, param ) &&
+ dist2 < distance2 )
+ {
+ distance2 = dist2;
+ u = param;
+ ok = true;
+ }
+ return ok;
+ }
+ else
+ {
+ for (int i = 0; i < 8; i++)
+ if (((CurveProjector*) myChildren[i])->IsOnCurve( point, dist2, param ) &&
+ dist2 < distance2 )
+ {
+ distance2 = dist2;
+ u = param;
+ ok = true;
+ }
+ }
+ return ok;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Initialize
+ */
+ //================================================================================
+
+ void CurveProjector::CurveSegment::Init(const gp_Pnt& pf,
+ const gp_Pnt& pl,
+ const double uf,
+ const double ul,
+ const double tol,
+ Handle(Geom_Curve)& curve )
+ {
+ _pFirst = pf;
+ _pLast = pl;
+ _curve = curve;
+ _length2 = pf.SquareDistance( pl );
+ _line.SetLocation( pf );
+ _line.SetDirection( gp_Vec( pf, pl ));
+ _chord2 = Max( _line. SquareDistance( curve->Value( uf + 0.25 * ( ul - uf ))),
+ Max( _line.SquareDistance( curve->Value( uf + 0.5 * ( ul - uf ))),
+ _line.SquareDistance( curve->Value( uf + 0.75 * ( ul - uf )))));
+ _chord2 *= ( 1.05 * 1.05 ); // +5%
+ _chord2 = Max( tol, _chord2 );
+ _chord = Sqrt( _chord2 );
+
+ Bnd_Box bb;
+ BndLib_Add3dCurve::Add( GeomAdaptor_Curve( curve, uf, ul ), tol, bb );
+ Add( bb.CornerMin() );
+ Add( bb.CornerMax() );
+ }
+
+ //================================================================================
+ /*!
+ * \brief Return true if a point is close to the curve segment
+ * \param [in] point - the point
+ * \param [out] distance2 - distance to the curve
+ * \param [out] u - parameter on the curve
+ * \return bool - is the point is close to the curve segment
+ */
+ //================================================================================
+
+ bool CurveProjector::CurveSegment::IsOn( const gp_XYZ& point, double & distance2, double & u )
+ {
+ distance2 = _line.SquareDistance( point );
+ if ( distance2 > _chord2 )
+ return false;
+
+ // check if the point projection falls into the segment range
+ {
+ gp_Vec edge( _pFirst, _pLast );
+ gp_Vec n1p ( _pFirst, point );
+ u = ( edge * n1p ) / _length2; // param [0,1] on the edge
+ if ( u < 0. )
+ {
+ if ( _pFirst.SquareDistance( point ) > _chord2 )
+ return false;
+ }
+ else if ( u > 1. )
+ {
+ if ( _pLast.SquareDistance( point ) > _chord2 )
+ return false;
+ }
+ }
+ gp_Pnt proj;
+ distance2 = ShapeAnalysis_Curve().Project( _curve, point, Precision::Confusion(),
+ proj, u, false );
+ distance2 *= distance2;
+ return true;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Check if the segment is in contact with a box
+ */
+ //================================================================================
+
+ bool CurveProjector::CurveSegment::IsInContact( const Bnd_B3d& bb )
+ {
+ if ( bb.IsOut( _line.Position(), /*isRay=*/true, _chord ))
+ return false;
+
+ gp_Ax1 axRev = _line.Position().Reversed();
+ axRev.SetLocation( _pLast );
+ return !bb.IsOut( axRev, /*isRay=*/true, _chord );
+ }
+
+ //================================================================================
+ //================================================================================
+