src/ADAPTFrontTrack/FrontTrack_Projector.cxx

   1 // Copyright (C) 2017-2020  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // File      : FrontTrack_Projector.cxx
  20 // Created   : Wed Apr 26 20:33:55 2017
  21 // Author    : Edward AGAPOV (eap)
  22
  23 #include "FrontTrack_Projector.hxx"
  24
  25 #include <BRepAdaptor_Curve.hxx>
  26 #include <BRepBndLib.hxx>
  27 #include <BRepTopAdaptor_FClass2d.hxx>
  28 #include <BRep_Tool.hxx>
  29 #include <ElCLib.hxx>
  30 #include <ElSLib.hxx>
  31 #include <GCPnts_UniformDeflection.hxx>
  32 #include <GeomAdaptor_Curve.hxx>
  33 #include <GeomLib_IsPlanarSurface.hxx>
  34 #include <ShapeAnalysis_Curve.hxx>
  35 #include <ShapeAnalysis_Surface.hxx>
  36 #include <TopExp.hxx>
  37 #include <TopoDS.hxx>
  38 #include <TopoDS_Edge.hxx>
  39 #include <TopoDS_Face.hxx>
  40 #include <TopoDS_Vertex.hxx>
  41 #include <gp_Circ.hxx>
  42 #include <gp_Cylinder.hxx>
  43 #include <gp_Dir.hxx>
  44 #include <gp_Pln.hxx>
  45 #include <gp_Pnt.hxx>
  46 #include <gp_Sphere.hxx>
  47 #include <gp_Vec.hxx>
  48
  49 #include <limits>
  50
  51 //-----------------------------------------------------------------------------
  52 /*!
  53  * \brief Root class of a projector of a point to a boundary shape
  54  */
  55 struct FT_RealProjector
  56 {
  57   virtual ~FT_RealProjector() {}
  58
  59   /*!
  60    * \brief Project a point to a boundary shape
  61    *  \param [in] point - the point to project
  62    *  \param [out] newSolution - position on the shape (U or UV) found during the projection
  63    *  \param [in] prevSolution - position already found during the projection of a neighbor point
  64    *  \return gp_Pnt - the projection point
  65    */
  66   virtual gp_Pnt project( const gp_Pnt& point,
  67                           double*       newSolution,
  68                           const double* prevSolution = 0) = 0;
  69
  70   /*!
  71    * \brief Project a point to a boundary shape and check if the projection is within
  72    *        the shape boundary
  73    *  \param [in] point - the point to project
  74    *  \param [in] maxDist2 - the maximal allowed square distance between point and projection
  75    *  \param [out] projection - the projection point
  76    *  \param [out] newSolution - position on the shape (U or UV) found during the projection
  77    *  \param [in] prevSolution - position already found during the projection of a neighbor point
  78    *  \return bool - false if the projection point lies out of the shape boundary or
  79    the distance the point and the projection is more than sqrt(maxDist2)
  80   */
  81   virtual bool projectAndClassify( const gp_Pnt& point,
  82                                    const double  maxDist2,
  83                                    gp_Pnt&       projection,
  84                                    double*       newSolution,
  85                                    const double* prevSolution = 0) = 0;
  86
  87   // return true if a previously found solution can be used to speed up the projection
  88
  89   virtual bool canUsePrevSolution() const { return false; }
  90
  91
  92   double _dist; // distance between the point being projected and its projection
  93 };
  94
  95 namespace // actual projection algorithms
  96 {
  97   const double theEPS = 1e-12;
  98
  99   //================================================================================
 100   /*!
 101    * \brief Projector to any curve
 102    */
 103   //================================================================================
 104
 105   struct CurveProjector : public FT_RealProjector
 106   {
 107     BRepAdaptor_Curve   _curve;
 108     double              _tol;
 109     ShapeAnalysis_Curve _projector;
 110     double              _uRange[2];
 111
 112     //-----------------------------------------------------------------------------
 113     CurveProjector( const TopoDS_Edge& e, const double tol ):
 114       _curve( e ), _tol( tol )
 115     {
 116       BRep_Tool::Range( e, _uRange[0], _uRange[1] );
 117     }
 118
 119     //-----------------------------------------------------------------------------
 120     // project a point to the curve
 121     virtual gp_Pnt project( const gp_Pnt& P,
 122                             double*       newSolution,
 123                             const double* prevSolution = 0)
 124     {
 125 #ifdef _DEBUG_
 126     std::cout << ".. project a point to the curve prevSolution = " << prevSolution << std::endl;
 127 #endif
 128       gp_Pnt         proj;
 129       Standard_Real param;
 130
 131       if ( prevSolution )
 132       {
 133         _dist = _projector.NextProject( prevSolution[0], _curve, P, _tol, proj, param );
 134       }
 135       else
 136       {
 137         _dist = _projector.Project( _curve, P, _tol, proj, param, false );
 138       }
 139 #ifdef _DEBUG_
 140     std::cout << "..    _dist : " << _dist << std::endl;
 141 #endif
 142       proj = _curve.Value( param );
 143
 144       newSolution[0] = param;
 145
 146       return proj;
 147     }
 148
 149     //-----------------------------------------------------------------------------
 150     // project a point to a curve and check if the projection is within the curve boundary
 151     virtual bool projectAndClassify( const gp_Pnt& point,
 152                                      const double  maxDist2,
 153                                      gp_Pnt&       projection,
 154                                      double*       newSolution,
 155                                      const double* prevSolution = 0)
 156     {
 157 #ifdef _DEBUG_
 158     std::cout << ".. project a point to a curve and check " << std::endl;
 159 #endif
 160       projection = project( point, newSolution, prevSolution );
 161       return ( _uRange[0] < newSolution[0] && newSolution[0] < _uRange[1] &&
 162                _dist * _dist < maxDist2 );
 163     }
 164
 165     //-----------------------------------------------------------------------------
 166     // return true if a previously found solution can be used to speed up the projection
 167     virtual bool canUsePrevSolution() const { return true; }
 168   };
 169
 170   //================================================================================
 171   /*!
 172    * \brief Projector to a straight curve. Don't project, classify only
 173    */
 174   //================================================================================
 175
 176   struct LineProjector : public FT_RealProjector
 177   {
 178     gp_Pnt _p0, _p1;
 179
 180     //-----------------------------------------------------------------------------
 181     LineProjector( TopoDS_Edge e )
 182     {
 183       e.Orientation( TopAbs_FORWARD );
 184       _p0 = BRep_Tool::Pnt( TopExp::FirstVertex( e ));
 185       _p1 = BRep_Tool::Pnt( TopExp::LastVertex ( e ));
 186     }
 187
 188     //-----------------------------------------------------------------------------
 189     // does nothing
 190     virtual gp_Pnt project( const gp_Pnt& P,
 191                             double*       newSolution,
 192                             const double* prevSolution = 0)
 193     {
 194       return P;
 195     }
 196     //-----------------------------------------------------------------------------
 197     // check if a point lies within the line segment
 198     virtual bool projectAndClassify( const gp_Pnt& point,
 199                                      const double  maxDist2,
 200                                      gp_Pnt&       projection,
 201                                      double*       newSolution,
 202                                      const double* prevSolution = 0)
 203     {
 204       gp_Vec edge( _p0, _p1 );
 205       gp_Vec p0p ( _p0, point  );
 206       double u = ( edge * p0p ) / edge.SquareMagnitude();  // param [0,1] on the edge
 207       projection = ( 1. - u ) * _p0.XYZ() + u * _p1.XYZ(); // projection of the point on the edge
 208       if ( u < 0 || 1 < u )
 209         return false;
 210
 211       // check distance
 212       return point.SquareDistance( projection ) < theEPS * theEPS;
 213     }
 214   };
 215
 216   //================================================================================
 217   /*!
 218    * \brief Projector to a circular edge
 219    */
 220   //================================================================================
 221
 222   struct CircleProjector : public FT_RealProjector
 223   {
 224     gp_Circ _circle;
 225     double _uRange[2];
 226
 227     //-----------------------------------------------------------------------------
 228     CircleProjector( const gp_Circ& c, const double f, const double l ):
 229       _circle( c )
 230     {
 231       _uRange[0] = f;
 232       _uRange[1] = l;
 233     }
 234
 235     //-----------------------------------------------------------------------------
 236     // project a point to the circle
 237     virtual gp_Pnt project( const gp_Pnt& P,
 238                             double*       newSolution,
 239                             const double* prevSolution = 0)
 240     {
 241       // assume that P is already on the the plane of circle, since
 242       // it is in the middle of two points lying on the circle
 243
 244       // move P to the circle
 245       const gp_Pnt& O = _circle.Location();
 246       gp_Vec radiusVec( O, P );
 247       double radius = radiusVec.Magnitude();
 248       if ( radius < std::numeric_limits<double>::min() )
 249         return P; // P in on the axe
 250
 251       gp_Pnt proj = O.Translated( radiusVec.Multiplied( _circle.Radius() / radius ));
 252
 253       _dist = _circle.Radius() - radius;
 254
 255       return proj;
 256     }
 257
 258     //-----------------------------------------------------------------------------
 259     // project and check if a projection lies within the circular edge
 260     virtual bool projectAndClassify( const gp_Pnt& point,
 261                                      const double  maxDist2,
 262                                      gp_Pnt&       projection,
 263                                      double*       newSolution,
 264                                      const double* prevSolution = 0)
 265     {
 266       _dist = -1;
 267       projection = project( point, newSolution );
 268       if ( _dist < 0 || // ?
 269            _dist * _dist > maxDist2 )
 270         return false;
 271
 272       newSolution[0] = ElCLib::Parameter( _circle, projection );
 273       return ( _uRange[0] < newSolution[0] && newSolution[0] < _uRange[1] );
 274     }
 275   };
 276
 277   //================================================================================
 278   /*!
 279    * \brief Projector to any surface
 280    */
 281   //================================================================================
 282
 283   struct SurfaceProjector : public FT_RealProjector
 284   {
 285     ShapeAnalysis_Surface    _projector;
 286     double                   _tol;
 287     BRepTopAdaptor_FClass2d* _classifier;
 288
 289     //-----------------------------------------------------------------------------
 290     SurfaceProjector( const TopoDS_Face& face, const double tol, BRepTopAdaptor_FClass2d* cls ):
 291       _projector( BRep_Tool::Surface( face )),
 292       _tol( tol ),
 293       _classifier( cls )
 294     {
 295     }
 296     //-----------------------------------------------------------------------------
 297     // delete _classifier
 298     ~SurfaceProjector()
 299     {
 300       delete _classifier;
 301     }
 302
 303     //-----------------------------------------------------------------------------
 304     // project a point to a surface
 305     virtual gp_Pnt project( const gp_Pnt& P,
 306                             double*       newSolution,
 307                             const double* prevSolution = 0)
 308     {
 309       gp_Pnt2d uv;
 310
 311       if ( prevSolution )
 312       {
 313         gp_Pnt2d prevUV( prevSolution[0], prevSolution[1] );
 314         uv = _projector.NextValueOfUV( prevUV, P, _tol );
 315       }
 316       else
 317       {
 318         uv = _projector.ValueOfUV( P, _tol );
 319       }
 320
 321       uv.Coord( newSolution[0], newSolution[1] );
 322
 323       gp_Pnt proj = _projector.Value( uv );
 324
 325       _dist = _projector.Gap();
 326
 327       return proj;
 328     }
 329
 330     //-----------------------------------------------------------------------------
 331     // project a point to a surface and check if the projection is within the surface boundary
 332     virtual bool projectAndClassify( const gp_Pnt& point,
 333                                      const double  maxDist2,
 334                                      gp_Pnt&       projection,
 335                                      double*       newSolution,
 336                                      const double* prevSolution = 0)
 337     {
 338       projection = project( point, newSolution, prevSolution );
 339       return ( _dist * _dist < maxDist2 )  &&  classify( newSolution );
 340     }
 341
 342     //-----------------------------------------------------------------------------
 343     // check if the projection is within the shape boundary
 344     bool classify( const double* newSolution )
 345     {
 346       TopAbs_State state = _classifier->Perform( gp_Pnt2d( newSolution[0], newSolution[1]) );
 347       return ( state != TopAbs_OUT );
 348     }
 349
 350     //-----------------------------------------------------------------------------
 351     // return true if a previously found solution can be used to speed up the projection
 352     virtual bool canUsePrevSolution() const { return true; }
 353   };
 354
 355   //================================================================================
 356   /*!
 357    * \brief Projector to a plane. Don't project, classify only
 358    */
 359   //================================================================================
 360
 361   struct PlaneProjector : public SurfaceProjector
 362   {
 363     gp_Pln _plane;
 364     bool   _isRealPlane; // false means that a surface is planar but parametrization is different
 365
 366     //-----------------------------------------------------------------------------
 367     PlaneProjector( const gp_Pln&            pln,
 368                     const TopoDS_Face&       face,
 369                     BRepTopAdaptor_FClass2d* cls,
 370                     bool                     isRealPlane=true):
 371       SurfaceProjector( face, 0, cls ),
 372       _plane( pln ),
 373       _isRealPlane( isRealPlane )
 374     {}
 375
 376     //-----------------------------------------------------------------------------
 377     // does nothing
 378     virtual gp_Pnt project( const gp_Pnt& P,
 379                             double*       newSolution,
 380                             const double* prevSolution = 0)
 381     {
 382       return P;
 383     }
 384     //-----------------------------------------------------------------------------
 385     // check if a point lies within the boundry of the planar face
 386     virtual bool projectAndClassify( const gp_Pnt& point,
 387                                      const double  maxDist2,
 388                                      gp_Pnt&       projection,
 389                                      double*       newSolution,
 390                                      const double* prevSolution = 0)
 391     {
 392       if ( _isRealPlane )
 393       {
 394         ElSLib::PlaneParameters( _plane.Position(), point, newSolution[0], newSolution[1]);
 395         projection = ElSLib::PlaneValue ( newSolution[0], newSolution[1], _plane.Position() );
 396         if ( projection.SquareDistance( point ) > theEPS * theEPS )
 397           return false;
 398
 399         return SurfaceProjector::classify( newSolution );
 400       }
 401       else
 402       {
 403         return SurfaceProjector::projectAndClassify( point, maxDist2, projection,
 404                                                      newSolution, prevSolution );
 405       }
 406     }
 407     //-----------------------------------------------------------------------------
 408     // return true if a previously found solution can be used to speed up the projection
 409     virtual bool canUsePrevSolution() const { return false; }
 410   };
 411
 412   //================================================================================
 413   /*!
 414    * \brief Projector to a cylinder
 415    */
 416   //================================================================================
 417
 418   struct CylinderProjector : public SurfaceProjector
 419   {
 420     gp_Cylinder _cylinder;
 421
 422     //-----------------------------------------------------------------------------
 423     CylinderProjector( const gp_Cylinder&       c,
 424                        const TopoDS_Face&       face,
 425                        BRepTopAdaptor_FClass2d* cls ):
 426       SurfaceProjector( face, 0, cls ),
 427       _cylinder( c )
 428     {}
 429
 430     //-----------------------------------------------------------------------------
 431     // project a point to the cylinder
 432     virtual gp_Pnt project( const gp_Pnt& P,
 433                             double*       newSolution,
 434                             const double* prevSolution = 0)
 435     {
 436       // project the point P to the cylinder axis -> Pp
 437       const gp_Pnt& O   = _cylinder.Position().Location();
 438       const gp_Dir& axe = _cylinder.Position().Direction();
 439       gp_Vec       trsl = gp_Vec( axe ).Multiplied( gp_Vec( O, P ).Dot( axe ));
 440       gp_Pnt       Pp   = O.Translated( trsl );
 441
 442       // move Pp to the cylinder
 443       gp_Vec radiusVec( Pp, P );
 444       double radius = radiusVec.Magnitude();
 445       if ( radius < std::numeric_limits<double>::min() )
 446         return P; // P in on the axe
 447
 448       gp_Pnt proj = Pp.Translated( radiusVec.Multiplied( _cylinder.Radius() / radius ));
 449
 450       _dist = _cylinder.Radius() - radius;
 451
 452       return proj;
 453     }
 454     //-----------------------------------------------------------------------------
 455     // project a point to the cylinder and check if the projection is within the surface boundary
 456     virtual bool projectAndClassify( const gp_Pnt& point,
 457                                      const double  maxDist2,
 458                                      gp_Pnt&       projection,
 459                                      double*       newSolution,
 460                                      const double* prevSolution = 0)
 461     {
 462       ElSLib::CylinderParameters( _cylinder.Position(), _cylinder.Radius(), point,
 463                                   newSolution[0], newSolution[1]);
 464       projection = ElSLib::CylinderValue( newSolution[0], newSolution[1],
 465                                           _cylinder.Position(), _cylinder.Radius() );
 466
 467       return ( _dist * _dist < maxDist2 )  &&  SurfaceProjector::classify( newSolution );
 468     }
 469     //-----------------------------------------------------------------------------
 470     // return true if a previously found solution can be used to speed up the projection
 471     virtual bool canUsePrevSolution() const { return false; }
 472   };
 473
 474   //================================================================================
 475   /*!
 476    * \brief Projector to a cone
 477    */
 478   //================================================================================
 479
 480   struct ConeProjector : public SurfaceProjector
 481   {
 482     gp_Cone _cone;
 483
 484     //-----------------------------------------------------------------------------
 485     ConeProjector( const gp_Cone&           c,
 486                    const TopoDS_Face&       face,
 487                    BRepTopAdaptor_FClass2d* cls ):
 488       SurfaceProjector( face, 0, cls ),
 489       _cone( c )
 490     {}
 491
 492     //-----------------------------------------------------------------------------
 493     // project a point to the cone
 494     virtual gp_Pnt project( const gp_Pnt& point,
 495                             double*       newSolution,
 496                             const double* prevSolution = 0)
 497     {
 498       ElSLib::ConeParameters( _cone.Position(), _cone.RefRadius(), _cone.SemiAngle(),
 499                               point, newSolution[0], newSolution[1]);
 500       gp_Pnt proj = ElSLib::ConeValue( newSolution[0], newSolution[1],
 501                                        _cone.Position(), _cone.RefRadius(), _cone.SemiAngle() );
 502       _dist = point.Distance( proj );
 503
 504       return proj;
 505     }
 506
 507     //-----------------------------------------------------------------------------
 508     // project a point to the cone and check if the projection is within the surface boundary
 509     virtual bool projectAndClassify( const gp_Pnt& point,
 510                                      const double  maxDist2,
 511                                      gp_Pnt&       projection,
 512                                      double*       newSolution,
 513                                      const double* prevSolution = 0)
 514     {
 515       projection = project( point, newSolution, prevSolution );
 516
 517       return ( _dist * _dist < maxDist2 )  &&  SurfaceProjector::classify( newSolution );
 518     }
 519     //-----------------------------------------------------------------------------
 520     // return true if a previously found solution can be used to speed up the projection
 521     virtual bool canUsePrevSolution() const { return false; }
 522   };
 523
 524   //================================================================================
 525   /*!
 526    * \brief Projector to a sphere
 527    */
 528   //================================================================================
 529
 530   struct SphereProjector : public SurfaceProjector
 531   {
 532     gp_Sphere _sphere;
 533
 534     //-----------------------------------------------------------------------------
 535     SphereProjector( const gp_Sphere&         s,
 536                      const TopoDS_Face&       face,
 537                      BRepTopAdaptor_FClass2d* cls ):
 538       SurfaceProjector( face, 0, cls ),
 539       _sphere( s )
 540     {}
 541
 542     //-----------------------------------------------------------------------------
 543     // project a point to the sphere
 544     virtual gp_Pnt project( const gp_Pnt& P,
 545                             double*       newSolution,
 546                             const double* prevSolution = 0)
 547     {
 548       // move Pp to the Sphere
 549       const gp_Pnt& O = _sphere.Location();
 550       gp_Vec radiusVec( O, P );
 551       double radius = radiusVec.Magnitude();
 552       if ( radius < std::numeric_limits<double>::min() )
 553         return P; // P is on O
 554
 555       gp_Pnt proj = O.Translated( radiusVec.Multiplied( _sphere.Radius() / radius ));
 556
 557       _dist = _sphere.Radius() - radius;
 558
 559       return proj;
 560     }
 561
 562     //-----------------------------------------------------------------------------
 563     // project a point to the sphere and check if the projection is within the surface boundary
 564     virtual bool projectAndClassify( const gp_Pnt& point,
 565                                      const double  maxDist2,
 566                                      gp_Pnt&       projection,
 567                                      double*       newSolution,
 568                                      const double* prevSolution = 0)
 569     {
 570       ElSLib::SphereParameters( _sphere.Position(), _sphere.Radius(), point,
 571                                   newSolution[0], newSolution[1]);
 572       projection = ElSLib::SphereValue( newSolution[0], newSolution[1],
 573                                         _sphere.Position(), _sphere.Radius() );
 574
 575       return ( _dist * _dist < maxDist2 )  &&  SurfaceProjector::classify( newSolution );
 576     }
 577     //-----------------------------------------------------------------------------
 578     // return true if a previously found solution can be used to speed up the projection
 579     virtual bool canUsePrevSolution() const { return false; }
 580   };
 581
 582   //================================================================================
 583   /*!
 584    * \brief Projector to a torus
 585    */
 586   //================================================================================
 587
 588   struct TorusProjector : public SurfaceProjector
 589   {
 590     gp_Torus _torus;
 591
 592     //-----------------------------------------------------------------------------
 593     TorusProjector( const gp_Torus&          t,
 594                     const TopoDS_Face&       face,
 595                     BRepTopAdaptor_FClass2d* cls ):
 596       SurfaceProjector( face, 0, cls ),
 597       _torus( t )
 598     {}
 599
 600     //-----------------------------------------------------------------------------
 601     // project a point to the torus
 602     virtual gp_Pnt project( const gp_Pnt& point,
 603                             double*       newSolution,
 604                             const double* prevSolution = 0)
 605     {
 606       ElSLib::TorusParameters( _torus.Position(), _torus.MajorRadius(), _torus.MinorRadius(),
 607                                point, newSolution[0], newSolution[1]);
 608       gp_Pnt proj = ElSLib::TorusValue( newSolution[0], newSolution[1],
 609                                         _torus.Position(), _torus.MajorRadius(), _torus.MinorRadius() );
 610       _dist = point.Distance( proj );
 611
 612       return proj;
 613     }
 614
 615     //-----------------------------------------------------------------------------
 616     // project a point to the torus and check if the projection is within the surface boundary
 617     virtual bool projectAndClassify( const gp_Pnt& point,
 618                                      const double  maxDist2,
 619                                      gp_Pnt&       projection,
 620                                      double*       newSolution,
 621                                      const double* prevSolution = 0)
 622     {
 623       projection = project( point, newSolution, prevSolution );
 624
 625       return ( _dist * _dist < maxDist2 )  &&  SurfaceProjector::classify( newSolution );
 626     }
 627     //-----------------------------------------------------------------------------
 628     // return true if a previously found solution can be used to speed up the projection
 629     virtual bool canUsePrevSolution() const { return false; }
 630   };
 631
 632   //================================================================================
 633   /*!
 634    * \brief Check if a curve can be considered straight
 635    */
 636   //================================================================================
 637
 638   bool isStraight( const GeomAdaptor_Curve& curve, const double tol )
 639   {
 640     // rough check: evaluate how far from a straight line connecting the curve ends
 641     // stand several internal points of the curve
 642
 643     const double  f = curve.FirstParameter();
 644     const double  l = curve.LastParameter();
 645     const gp_Pnt pf = curve.Value( f );
 646     const gp_Pnt pl = curve.Value( l );
 647     const gp_Vec lineVec( pf, pl );
 648     const double lineLen2 = lineVec.SquareMagnitude();
 649     if ( lineLen2 < std::numeric_limits< double >::min() )
 650       return false; // E seems closed
 651
 652     const double nbSamples = 7;
 653     for ( int i = 0; i < nbSamples; ++i )
 654     {
 655       const double  r = ( i + 1 ) / nbSamples;
 656       const gp_Pnt pi = curve.Value( f * r + l * ( 1 - r ));
 657       const gp_Vec vi( pf, pi );
 658       const double h2 = lineVec.Crossed( vi ).SquareMagnitude() / lineLen2;
 659       if ( h2 > tol * tol )
 660         return false;
 661     }
 662
 663     // thorough check
 664     GCPnts_UniformDeflection divider( curve, tol );
 665     return ( divider.IsDone() && divider.NbPoints() < 3 );
 666   }
 667 }
 668
 669 //================================================================================
 670 /*!
 671  * \brief Initialize with a boundary shape
 672  */
 673 //================================================================================
 674
 675 FT_Projector::FT_Projector(const TopoDS_Shape& shape)
 676 {
 677   _realProjector = 0;
 678   setBoundaryShape( shape );
 679   _tryWOPrevSolution = false;
 680 }
 681
 682 //================================================================================
 683 /*!
 684  * \brief Copy another projector
 685  */
 686 //================================================================================
 687
 688 FT_Projector::FT_Projector(const FT_Projector& other)
 689 {
 690   _realProjector = 0;
 691   _shape = other._shape;
 692   _bndBox = other._bndBox;
 693   _tryWOPrevSolution = false;
 694 }
 695
 696 //================================================================================
 697 /*!
 698  * \brief Destructor. Delete _realProjector
 699  */
 700 //================================================================================
 701
 702 FT_Projector::~FT_Projector()
 703 {
 704   delete _realProjector;
 705 }
 706
 707 //================================================================================
 708 /*!
 709  * \brief Initialize with a boundary shape. Compute the bounding box
 710  */
 711 //================================================================================
 712
 713 void FT_Projector::setBoundaryShape(const TopoDS_Shape& shape)
 714 {
 715   delete _realProjector; _realProjector = 0;
 716   _shape = shape;
 717   if ( shape.IsNull() )
 718     return;
 719
 720   BRepBndLib::Add( shape, _bndBox );
 721   _bndBox.Enlarge( 1e-5 * sqrt( _bndBox.SquareExtent() ));
 722 }
 723
 724 //================================================================================
 725 /*!
 726  * \brief Create a real projector
 727  */
 728 //================================================================================
 729
 730 void FT_Projector::prepareForProjection()
 731 {
 732   if ( _shape.IsNull() || _realProjector )
 733     return;
 734
 735   if ( _shape.ShapeType() == TopAbs_EDGE )
 736   {
 737     const TopoDS_Edge& edge = TopoDS::Edge( _shape );
 738
 739     double tol = 1e-6 * sqrt( _bndBox.SquareExtent() );
 740
 741     double f,l;
 742     Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, f,l );
 743     if ( curve.IsNull() )
 744       return; // degenerated edge
 745
 746     GeomAdaptor_Curve acurve( curve, f, l );
 747     switch ( acurve.GetType() )
 748     {
 749     case GeomAbs_Line:
 750       _realProjector = new LineProjector( edge );
 751       break;
 752     case GeomAbs_Circle:
 753       _realProjector = new CircleProjector( acurve.Circle(), f, l );
 754       break;
 755     case GeomAbs_BezierCurve:
 756     case GeomAbs_BSplineCurve:
 757     case GeomAbs_OffsetCurve:
 758     case GeomAbs_OtherCurve:
 759       if ( isStraight( acurve, tol ))
 760       {
 761         _realProjector = new LineProjector( edge );
 762         break;
 763       }
 764     case GeomAbs_Ellipse:
 765     case GeomAbs_Hyperbola:
 766     case GeomAbs_Parabola:
 767       _realProjector = new CurveProjector( edge, tol );
 768     }
 769   }
 770   else if ( _shape.ShapeType() == TopAbs_FACE )
 771   {
 772     TopoDS_Face face = TopoDS::Face( _shape );
 773
 774     Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
 775     if ( surface.IsNull() )
 776       return;
 777
 778     GeomAdaptor_Surface asurface( surface );
 779     Standard_Real tol   = BRep_Tool::Tolerance( face );
 780     Standard_Real toluv = Min( asurface.UResolution( tol ), asurface.VResolution( tol ));
 781     BRepTopAdaptor_FClass2d* classifier = new BRepTopAdaptor_FClass2d( face, toluv );
 782
 783     switch ( asurface.GetType() )
 784     {
 785     case GeomAbs_Plane:
 786       _realProjector = new PlaneProjector( asurface.Plane(), face, classifier );
 787       break;
 788     case GeomAbs_Cylinder:
 789       _realProjector = new CylinderProjector( asurface.Cylinder(), face, classifier );
 790       break;
 791     case GeomAbs_Sphere:
 792       _realProjector = new SphereProjector( asurface.Sphere(), face, classifier );
 793       break;
 794     case GeomAbs_Cone:
 795       _realProjector = new ConeProjector( asurface.Cone(), face, classifier );
 796       break;
 797     case GeomAbs_Torus:
 798       _realProjector = new TorusProjector( asurface.Torus(), face, classifier );
 799       break;
 800     case GeomAbs_BezierSurface:
 801     case GeomAbs_BSplineSurface:
 802     case GeomAbs_SurfaceOfRevolution:
 803     case GeomAbs_SurfaceOfExtrusion:
 804     case GeomAbs_OffsetSurface:
 805     case GeomAbs_OtherSurface:
 806       GeomLib_IsPlanarSurface isPlaneCheck( surface, tol );
 807       if ( isPlaneCheck.IsPlanar() )
 808       {
 809         _realProjector = new PlaneProjector( isPlaneCheck.Plan(), face, classifier,
 810                                              /*isRealPlane=*/false);
 811       }
 812       else
 813       {
 814         _realProjector = new SurfaceProjector( face, tol, classifier );
 815       }
 816       break;
 817     }
 818
 819     if ( !_realProjector )
 820       delete classifier;
 821   }
 822 }
 823
 824 //================================================================================
 825 /*!
 826  * \brief Return true if projection is not needed
 827  */
 828 //================================================================================
 829
 830 bool FT_Projector::isPlanarBoundary() const
 831 {
 832   return ( dynamic_cast< LineProjector*  >( _realProjector ) ||
 833            dynamic_cast< PlaneProjector* >( _realProjector ) );
 834 }
 835
 836 //================================================================================
 837 /*!
 838  * \brief Check if a point lies on the boundary shape
 839  *  \param [in] point - the point to check
 840  *  \param [in] tol2 - a square tolerance allowing to decide whether a point is on the shape
 841  *  \param [in] newSolution - position on the shape (U or UV) of the point found
 842  *         during projecting
 843  *  \param [in] prevSolution - position on the shape (U or UV) of a neighbor point
 844  *  \return bool - \c true if the point lies on the boundary shape
 845  *
 846  * This method is used to select a shape by checking if all neighbor nodes of a node to move
 847  * lie on a shape.
 848  */
 849 //================================================================================
 850
 851 bool FT_Projector::isOnShape( const gp_Pnt& point,
 852                               const double  tol2,
 853                               double*       newSolution,
 854                               const double* prevSolution)
 855 {
 856   if ( _bndBox.IsOut( point ) || !_realProjector )
 857     return false;
 858
 859   gp_Pnt proj;
 860   if ( isPlanarBoundary() )
 861     return projectAndClassify( point, tol2, proj, newSolution, prevSolution );
 862
 863   return project( point, tol2, proj, newSolution, prevSolution );
 864 }
 865
 866 //================================================================================
 867 /*!
 868  * \brief Project a point to the boundary shape
 869  *  \param [in] point - the point to project
 870  *  \param [in] maxDist2 - the maximal square distance between the point and the projection
 871  *  \param [out] projection - the projection
 872  *  \param [out] newSolution - position on the shape (U or UV) of the point found
 873  *         during projecting
 874  *  \param [in] prevSolution - already found position on the shape (U or UV) of a neighbor point
 875  *  \return bool - false if the distance between the point and the projection
 876  *         is more than sqrt(maxDist2)
 877  *
 878  * This method is used to project a node in the case where only one shape is found by name
 879  */
 880 //================================================================================
 881
 882 bool FT_Projector::project( const gp_Pnt& point,
 883                             const double  maxDist2,
 884                             gp_Pnt&       projection,
 885                             double*       newSolution,
 886                             const double* prevSolution)
 887 {
 888   if ( !_realProjector )
 889     return false;
 890
 891   _realProjector->_dist = 1e100;
 892   projection = _realProjector->project( point, newSolution, prevSolution );
 893
 894   bool ok = ( _realProjector->_dist * _realProjector->_dist < maxDist2 );
 895   if ( !ok && _tryWOPrevSolution && prevSolution )
 896   {
 897     projection = _realProjector->project( point, newSolution );
 898     ok = ( _realProjector->_dist * _realProjector->_dist < maxDist2 );
 899   }
 900   return ok;
 901 }
 902
 903 //================================================================================
 904 /*!
 905  * \brief Project a point to the boundary shape and check if the projection lies within
 906  *        the shape boundary
 907  *  \param [in] point - the point to project
 908  *  \param [in] maxDist2 - the maximal square distance between the point and the projection
 909  *  \param [out] projection - the projection
 910  *  \param [out] newSolution - position on the shape (U or UV) of the point found
 911  *         during projecting
 912  *  \param [in] prevSolution - already found position on the shape (U or UV) of a neighbor point
 913  *  \return bool - false if the projection point lies out of the shape boundary or
 914  *          the distance between the point and the projection is more than sqrt(maxDist2)
 915  *
 916  * This method is used to project a node in the case where several shapes are selected for
 917  * projection of a node group
 918  */
 919 //================================================================================
 920
 921 bool FT_Projector::projectAndClassify( const gp_Pnt& point,
 922                                        const double  maxDist2,
 923                                        gp_Pnt&       projection,
 924                                        double*       newSolution,
 925                                        const double* prevSolution)
 926 {
 927   if ( _bndBox.IsOut( point ) || !_realProjector )
 928     return false;
 929
 930   bool ok = _realProjector->projectAndClassify( point, maxDist2, projection,
 931                                                 newSolution, prevSolution );
 932   if ( !ok && _tryWOPrevSolution && prevSolution )
 933     ok = _realProjector->projectAndClassify( point, maxDist2, projection, newSolution );
 934
 935   return ok;
 936 }
 937
 938 //================================================================================
 939 /*!
 940  * \brief Return true if a previously found solution can be used to speed up the projection
 941  */
 942 //================================================================================
 943
 944 bool FT_Projector::canUsePrevSolution() const
 945 {
 946   return ( _realProjector && _realProjector->canUsePrevSolution() );
 947 }