Salome HOME
Copyright update 2022
[modules/homard.git] / src / FrontTrack / FrontTrack_Projector.cxx
1 // Copyright (C) 2017-2022  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       _dist = point.Distance( projection );
467
468       return ( _dist * _dist < maxDist2 )  &&  SurfaceProjector::classify( newSolution );
469     }
470     //-----------------------------------------------------------------------------
471     // return true if a previously found solution can be used to speed up the projection
472     virtual bool canUsePrevSolution() const { return false; }
473   };
474
475   //================================================================================
476   /*!
477    * \brief Projector to a cone
478    */
479   //================================================================================
480
481   struct ConeProjector : public SurfaceProjector
482   {
483     gp_Cone _cone;
484
485     //-----------------------------------------------------------------------------
486     ConeProjector( const gp_Cone&           c,
487                    const TopoDS_Face&       face,
488                    BRepTopAdaptor_FClass2d* cls ):
489       SurfaceProjector( face, 0, cls ),
490       _cone( c )
491     {}
492
493     //-----------------------------------------------------------------------------
494     // project a point to the cone
495     virtual gp_Pnt project( const gp_Pnt& point,
496                             double*       newSolution,
497                             const double* prevSolution = 0)
498     {
499       ElSLib::ConeParameters( _cone.Position(), _cone.RefRadius(), _cone.SemiAngle(),
500                               point, newSolution[0], newSolution[1]);
501       gp_Pnt proj = ElSLib::ConeValue( newSolution[0], newSolution[1],
502                                        _cone.Position(), _cone.RefRadius(), _cone.SemiAngle() );
503       _dist = point.Distance( proj );
504
505       return proj;
506     }
507
508     //-----------------------------------------------------------------------------
509     // project a point to the cone and check if the projection is within the surface boundary
510     virtual bool projectAndClassify( const gp_Pnt& point,
511                                      const double  maxDist2,
512                                      gp_Pnt&       projection,
513                                      double*       newSolution,
514                                      const double* prevSolution = 0)
515     {
516       projection = project( point, newSolution, prevSolution );
517
518       return ( _dist * _dist < maxDist2 )  &&  SurfaceProjector::classify( newSolution );
519     }
520     //-----------------------------------------------------------------------------
521     // return true if a previously found solution can be used to speed up the projection
522     virtual bool canUsePrevSolution() const { return false; }
523   };
524
525   //================================================================================
526   /*!
527    * \brief Projector to a sphere
528    */
529   //================================================================================
530
531   struct SphereProjector : public SurfaceProjector
532   {
533     gp_Sphere _sphere;
534
535     //-----------------------------------------------------------------------------
536     SphereProjector( const gp_Sphere&         s,
537                      const TopoDS_Face&       face,
538                      BRepTopAdaptor_FClass2d* cls ):
539       SurfaceProjector( face, 0, cls ),
540       _sphere( s )
541     {}
542
543     //-----------------------------------------------------------------------------
544     // project a point to the sphere
545     virtual gp_Pnt project( const gp_Pnt& P,
546                             double*       newSolution,
547                             const double* prevSolution = 0)
548     {
549       // move Pp to the Sphere
550       const gp_Pnt& O = _sphere.Location();
551       gp_Vec radiusVec( O, P );
552       double radius = radiusVec.Magnitude();
553       if ( radius < std::numeric_limits<double>::min() )
554         return P; // P is on O
555
556       gp_Pnt proj = O.Translated( radiusVec.Multiplied( _sphere.Radius() / radius ));
557
558       _dist = _sphere.Radius() - radius;
559
560       return proj;
561     }
562
563     //-----------------------------------------------------------------------------
564     // project a point to the sphere and check if the projection is within the surface boundary
565     virtual bool projectAndClassify( const gp_Pnt& point,
566                                      const double  maxDist2,
567                                      gp_Pnt&       projection,
568                                      double*       newSolution,
569                                      const double* prevSolution = 0)
570     {
571       ElSLib::SphereParameters( _sphere.Position(), _sphere.Radius(), point,
572                                   newSolution[0], newSolution[1]);
573       projection = ElSLib::SphereValue( newSolution[0], newSolution[1],
574                                         _sphere.Position(), _sphere.Radius() );
575
576       return ( _dist * _dist < maxDist2 )  &&  SurfaceProjector::classify( newSolution );
577     }
578     //-----------------------------------------------------------------------------
579     // return true if a previously found solution can be used to speed up the projection
580     virtual bool canUsePrevSolution() const { return false; }
581   };
582
583   //================================================================================
584   /*!
585    * \brief Projector to a torus
586    */
587   //================================================================================
588
589   struct TorusProjector : public SurfaceProjector
590   {
591     gp_Torus _torus;
592
593     //-----------------------------------------------------------------------------
594     TorusProjector( const gp_Torus&          t,
595                     const TopoDS_Face&       face,
596                     BRepTopAdaptor_FClass2d* cls ):
597       SurfaceProjector( face, 0, cls ),
598       _torus( t )
599     {}
600
601     //-----------------------------------------------------------------------------
602     // project a point to the torus
603     virtual gp_Pnt project( const gp_Pnt& point,
604                             double*       newSolution,
605                             const double* prevSolution = 0)
606     {
607       ElSLib::TorusParameters( _torus.Position(), _torus.MajorRadius(), _torus.MinorRadius(),
608                                point, newSolution[0], newSolution[1]);
609       gp_Pnt proj = ElSLib::TorusValue( newSolution[0], newSolution[1],
610                                         _torus.Position(), _torus.MajorRadius(), _torus.MinorRadius() );
611       _dist = point.Distance( proj );
612
613       return proj;
614     }
615
616     //-----------------------------------------------------------------------------
617     // project a point to the torus and check if the projection is within the surface boundary
618     virtual bool projectAndClassify( const gp_Pnt& point,
619                                      const double  maxDist2,
620                                      gp_Pnt&       projection,
621                                      double*       newSolution,
622                                      const double* prevSolution = 0)
623     {
624       projection = project( point, newSolution, prevSolution );
625
626       return ( _dist * _dist < maxDist2 )  &&  SurfaceProjector::classify( newSolution );
627     }
628     //-----------------------------------------------------------------------------
629     // return true if a previously found solution can be used to speed up the projection
630     virtual bool canUsePrevSolution() const { return false; }
631   };
632
633   //================================================================================
634   /*!
635    * \brief Check if a curve can be considered straight
636    */
637   //================================================================================
638
639   bool isStraight( const GeomAdaptor_Curve& curve, const double tol )
640   {
641     // rough check: evaluate how far from a straight line connecting the curve ends
642     // stand several internal points of the curve
643
644     const double  f = curve.FirstParameter();
645     const double  l = curve.LastParameter();
646     const gp_Pnt pf = curve.Value( f );
647     const gp_Pnt pl = curve.Value( l );
648     const gp_Vec lineVec( pf, pl );
649     const double lineLen2 = lineVec.SquareMagnitude();
650     if ( lineLen2 < std::numeric_limits< double >::min() )
651       return false; // E seems closed
652
653     const double nbSamples = 7;
654     for ( int i = 0; i < nbSamples; ++i )
655     {
656       const double  r = ( i + 1 ) / nbSamples;
657       const gp_Pnt pi = curve.Value( f * r + l * ( 1 - r ));
658       const gp_Vec vi( pf, pi );
659       const double h2 = lineVec.Crossed( vi ).SquareMagnitude() / lineLen2;
660       if ( h2 > tol * tol )
661         return false;
662     }
663
664     // thorough check
665     GCPnts_UniformDeflection divider( curve, tol );
666     return ( divider.IsDone() && divider.NbPoints() < 3 );
667   }
668 }
669
670 //================================================================================
671 /*!
672  * \brief Initialize with a boundary shape
673  */
674 //================================================================================
675
676 FT_Projector::FT_Projector(const TopoDS_Shape& shape)
677 {
678   _realProjector = 0;
679   setBoundaryShape( shape );
680   _tryWOPrevSolution = false;
681 }
682
683 //================================================================================
684 /*!
685  * \brief Copy another projector
686  */
687 //================================================================================
688
689 FT_Projector::FT_Projector(const FT_Projector& other)
690 {
691   _realProjector = 0;
692   _shape = other._shape;
693   _bndBox = other._bndBox;
694   _tryWOPrevSolution = false;
695 }
696
697 //================================================================================
698 /*!
699  * \brief Destructor. Delete _realProjector
700  */
701 //================================================================================
702
703 FT_Projector::~FT_Projector()
704 {
705   delete _realProjector;
706 }
707
708 //================================================================================
709 /*!
710  * \brief Initialize with a boundary shape. Compute the bounding box
711  */
712 //================================================================================
713
714 void FT_Projector::setBoundaryShape(const TopoDS_Shape& shape)
715 {
716   delete _realProjector; _realProjector = 0;
717   _shape = shape;
718   if ( shape.IsNull() )
719     return;
720
721   BRepBndLib::Add( shape, _bndBox );
722   _bndBox.Enlarge( 1e-5 * sqrt( _bndBox.SquareExtent() ));
723 }
724
725 //================================================================================
726 /*!
727  * \brief Create a real projector
728  */
729 //================================================================================
730
731 void FT_Projector::prepareForProjection()
732 {
733   if ( _shape.IsNull() || _realProjector )
734     return;
735
736   if ( _shape.ShapeType() == TopAbs_EDGE )
737   {
738     const TopoDS_Edge& edge = TopoDS::Edge( _shape );
739
740     double tol = 1e-6 * sqrt( _bndBox.SquareExtent() );
741
742     double f,l;
743     Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, f,l );
744     if ( curve.IsNull() )
745       return; // degenerated edge
746
747     GeomAdaptor_Curve acurve( curve, f, l );
748     switch ( acurve.GetType() )
749     {
750     case GeomAbs_Line:
751       _realProjector = new LineProjector( edge );
752       break;
753     case GeomAbs_Circle:
754       _realProjector = new CircleProjector( acurve.Circle(), f, l );
755       break;
756     case GeomAbs_BezierCurve:
757     case GeomAbs_BSplineCurve:
758     case GeomAbs_OffsetCurve:
759     case GeomAbs_OtherCurve:
760       if ( isStraight( acurve, tol ))
761       {
762         _realProjector = new LineProjector( edge );
763         break;
764       }
765     case GeomAbs_Ellipse:
766     case GeomAbs_Hyperbola:
767     case GeomAbs_Parabola:
768       _realProjector = new CurveProjector( edge, tol );
769     }
770   }
771   else if ( _shape.ShapeType() == TopAbs_FACE )
772   {
773     TopoDS_Face face = TopoDS::Face( _shape );
774
775     Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
776     if ( surface.IsNull() )
777       return;
778
779     GeomAdaptor_Surface asurface( surface );
780     Standard_Real tol   = BRep_Tool::Tolerance( face );
781     Standard_Real toluv = Min( asurface.UResolution( tol ), asurface.VResolution( tol ));
782     BRepTopAdaptor_FClass2d* classifier = new BRepTopAdaptor_FClass2d( face, toluv );
783
784     switch ( asurface.GetType() )
785     {
786     case GeomAbs_Plane:
787       _realProjector = new PlaneProjector( asurface.Plane(), face, classifier );
788       break;
789     case GeomAbs_Cylinder:
790       _realProjector = new CylinderProjector( asurface.Cylinder(), face, classifier );
791       break;
792     case GeomAbs_Sphere:
793       _realProjector = new SphereProjector( asurface.Sphere(), face, classifier );
794       break;
795     case GeomAbs_Cone:
796       _realProjector = new ConeProjector( asurface.Cone(), face, classifier );
797       break;
798     case GeomAbs_Torus:
799       _realProjector = new TorusProjector( asurface.Torus(), face, classifier );
800       break;
801     case GeomAbs_BezierSurface:
802     case GeomAbs_BSplineSurface:
803     case GeomAbs_SurfaceOfRevolution:
804     case GeomAbs_SurfaceOfExtrusion:
805     case GeomAbs_OffsetSurface:
806     case GeomAbs_OtherSurface:
807       GeomLib_IsPlanarSurface isPlaneCheck( surface, tol );
808       if ( isPlaneCheck.IsPlanar() )
809       {
810         _realProjector = new PlaneProjector( isPlaneCheck.Plan(), face, classifier,
811                                              /*isRealPlane=*/false);
812       }
813       else
814       {
815         _realProjector = new SurfaceProjector( face, tol, classifier );
816       }
817       break;
818     }
819
820     if ( !_realProjector )
821       delete classifier;
822   }
823 }
824
825 //================================================================================
826 /*!
827  * \brief Return true if projection is not needed
828  */
829 //================================================================================
830
831 bool FT_Projector::isPlanarBoundary() const
832 {
833   return ( dynamic_cast< LineProjector*  >( _realProjector ) ||
834            dynamic_cast< PlaneProjector* >( _realProjector ) );
835 }
836
837 //================================================================================
838 /*!
839  * \brief Check if a point lies on the boundary shape
840  *  \param [in] point - the point to check
841  *  \param [in] tol2 - a square tolerance allowing to decide whether a point is on the shape
842  *  \param [in] newSolution - position on the shape (U or UV) of the point found
843  *         during projecting
844  *  \param [in] prevSolution - position on the shape (U or UV) of a neighbor point
845  *  \return bool - \c true if the point lies on the boundary shape
846  *
847  * This method is used to select a shape by checking if all neighbor nodes of a node to move
848  * lie on a shape.
849  */
850 //================================================================================
851
852 bool FT_Projector::isOnShape( const gp_Pnt& point,
853                               const double  tol2,
854                               double*       newSolution,
855                               const double* prevSolution)
856 {
857   if ( _bndBox.IsOut( point ) || !_realProjector )
858     return false;
859
860   gp_Pnt proj;
861   if ( isPlanarBoundary() )
862     return projectAndClassify( point, tol2, proj, newSolution, prevSolution );
863
864   return project( point, tol2, proj, newSolution, prevSolution );
865 }
866
867 //================================================================================
868 /*!
869  * \brief Project a point to the boundary shape
870  *  \param [in] point - the point to project
871  *  \param [in] maxDist2 - the maximal square distance between the point and the projection
872  *  \param [out] projection - the projection
873  *  \param [out] newSolution - position on the shape (U or UV) of the point found
874  *         during projecting
875  *  \param [in] prevSolution - already found position on the shape (U or UV) of a neighbor point
876  *  \return bool - false if the distance between the point and the projection
877  *         is more than sqrt(maxDist2)
878  *
879  * This method is used to project a node in the case where only one shape is found by name
880  */
881 //================================================================================
882
883 bool FT_Projector::project( const gp_Pnt& point,
884                             const double  maxDist2,
885                             gp_Pnt&       projection,
886                             double*       newSolution,
887                             const double* prevSolution)
888 {
889   if ( !_realProjector )
890     return false;
891
892   _realProjector->_dist = 1e100;
893   projection = _realProjector->project( point, newSolution, prevSolution );
894
895   bool ok = ( _realProjector->_dist * _realProjector->_dist < maxDist2 );
896   if ( !ok && _tryWOPrevSolution && prevSolution )
897   {
898     projection = _realProjector->project( point, newSolution );
899     ok = ( _realProjector->_dist * _realProjector->_dist < maxDist2 );
900   }
901   return ok;
902 }
903
904 //================================================================================
905 /*!
906  * \brief Project a point to the boundary shape and check if the projection lies within
907  *        the shape boundary
908  *  \param [in] point - the point to project
909  *  \param [in] maxDist2 - the maximal square distance between the point and the projection
910  *  \param [out] projection - the projection
911  *  \param [out] newSolution - position on the shape (U or UV) of the point found
912  *         during projecting
913  *  \param [in] prevSolution - already found position on the shape (U or UV) of a neighbor point
914  *  \return bool - false if the projection point lies out of the shape boundary or
915  *          the distance between the point and the projection is more than sqrt(maxDist2)
916  *
917  * This method is used to project a node in the case where several shapes are selected for
918  * projection of a node group
919  */
920 //================================================================================
921
922 bool FT_Projector::projectAndClassify( const gp_Pnt& point,
923                                        const double  maxDist2,
924                                        gp_Pnt&       projection,
925                                        double*       newSolution,
926                                        const double* prevSolution)
927 {
928   if ( _bndBox.IsOut( point ) || !_realProjector )
929     return false;
930
931   bool ok = _realProjector->projectAndClassify( point, maxDist2, projection,
932                                                 newSolution, prevSolution );
933   if ( !ok && _tryWOPrevSolution && prevSolution )
934     ok = _realProjector->projectAndClassify( point, maxDist2, projection, newSolution );
935
936   return ok;
937 }
938
939 //================================================================================
940 /*!
941  * \brief Return true if a previously found solution can be used to speed up the projection
942  */
943 //================================================================================
944
945 bool FT_Projector::canUsePrevSolution() const
946 {
947   return ( _realProjector && _realProjector->canUsePrevSolution() );
948 }