src/IntPatch/IntPatch_ALineToWLine.cxx

   1 // Created on: 1993-11-26
   2 // Created by: Modelistation
   3 // Copyright (c) 1993-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and / or modify it
   9 // under the terms of the GNU Lesser General Public version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <IntPatch_ALineToWLine.ixx>
  18
  19 #include <IntSurf_LineOn2S.hxx>
  20 #include <IntSurf_PntOn2S.hxx>
  21 #include <IntSurf_TypeTrans.hxx>
  22 #include <IntSurf_Situation.hxx>
  23
  24 #include <TColStd_Array1OfReal.hxx>
  25 #include <TColStd_Array1OfInteger.hxx>
  26
  27 #include <Precision.hxx>
  28
  29 #include <gp_Pnt.hxx>
  30 #include <gp_Vec.hxx>
  31
  32 #include <Adaptor2d_HCurve2d.hxx>
  33 #include <GeomAbs_SurfaceType.hxx>
  34 #include <gp_Pnt2d.hxx>
  35 #include <gp_Vec2d.hxx>
  36 #include <IntAna2d_AnaIntersection.hxx>
  37 #include <gp_Lin2d.hxx>
  38 #include <IntAna2d_IntPoint.hxx>
  39 #include <gp_Cone.hxx>
  40
  41
  42 static
  43   gp_Pnt DefineDU(const Handle(IntPatch_ALine)& aline,
  44                   const Standard_Real U,
  45                   Standard_Real& DU,
  46                   const Standard_Real CurvDef,
  47                   const Standard_Real AngDef);
  48
  49 static
  50   Standard_Boolean SameCurve(const Handle_Adaptor2d_HCurve2d& C1,
  51                              const Handle_Adaptor2d_HCurve2d& C2);
  52
  53 static
  54   void RecadreMemePeriode(Standard_Real& u1,
  55                           Standard_Real& u2,
  56                           const Standard_Real anu1,
  57                           const Standard_Real anu2);
  58
  59 static
  60   void CorrectFirstPartOfLine(Handle(IntSurf_LineOn2S)& LinOn2S,
  61                               const Standard_Real ref_u1,
  62                               const Standard_Real ref_u2,
  63                               Standard_Real& new_u1,
  64                               Standard_Real& new_u2);
  65
  66 //
  67 static
  68   void RefineParameters(const Handle(IntPatch_ALine)& aALine,
  69                       const Standard_Real aTb,
  70                       const Standard_Real aTe,
  71                       const Standard_Real aTx,
  72                       const Standard_Integer iDir,
  73                       const IntSurf_Quadric& aQuadric,
  74                       const Standard_Real aTol3D,
  75                       Standard_Real& aUx,
  76                       Standard_Real& aVx);
  77 static
  78   Standard_Boolean FindNearParameter(const Handle(IntPatch_ALine)& aALine,
  79                                      const Standard_Real aTx,
  80                                      const Standard_Integer iDir,
  81                                      const Standard_Real aTol3D,
  82                                      Standard_Real& aT1);
  83 static
  84   Standard_Boolean IsApex(const IntSurf_Quadric& aQuadric,
  85                           const Standard_Real aVx,
  86                           const Standard_Real aTol3D);
  87 //
  88
  89 //--  Je ne sais pas faire mieux et ca m'ennerve.  lbr.
  90 //=======================================================================
  91 //function : IntPatch_ALineToWLine
  92 //purpose  :
  93 //=======================================================================
  94   IntPatch_ALineToWLine::IntPatch_ALineToWLine(const IntSurf_Quadric& Quad1,
  95                                                const IntSurf_Quadric& Quad2)
  96 :
  97   quad1(Quad1),
  98   quad2(Quad2),
  99   deflectionmax(0.01),
 100   nbpointsmax(200),
 101   type(0),
 102   myTolParam(1.e-12),
 103   myTolOpenDomain(1.e-9),
 104   myTolTransition(1.e-8)
 105 {
 106   myTol3D=Precision::Confusion();
 107 }
 108 //=======================================================================
 109 //function : IntPatch_ALineToWLine
 110 //purpose  :
 111 //=======================================================================
 112   IntPatch_ALineToWLine::IntPatch_ALineToWLine(const IntSurf_Quadric& Quad1,
 113                                                const IntSurf_Quadric& Quad2,
 114                                                const Standard_Real    Deflection,
 115                                                const Standard_Real    ,
 116                                                const Standard_Integer NbMaxPoints)
 117 :
 118   quad1(Quad1),
 119   quad2(Quad2),
 120   deflectionmax(Deflection),
 121   nbpointsmax(NbMaxPoints),
 122   myTolParam(1.e-12),
 123   myTolOpenDomain(1.e-9),
 124   myTolTransition(1.e-8)
 125 {
 126   myTol3D=Precision::Confusion();
 127   type = 0;
 128 }
 129 //=======================================================================
 130 //function : SetTol3D
 131 //purpose  :
 132 //=======================================================================
 133   void IntPatch_ALineToWLine::SetTol3D(const Standard_Real aTol)
 134 {
 135   myTol3D = aTol;
 136 }
 137 //=======================================================================
 138 //function : Tol3D
 139 //purpose  :
 140 //=======================================================================
 141   Standard_Real IntPatch_ALineToWLine::Tol3D()const
 142 {
 143   return myTol3D;
 144 }
 145 //=======================================================================
 146 //function : SetTolTransition
 147 //purpose  :
 148 //=======================================================================
 149   void IntPatch_ALineToWLine::SetTolTransition(const Standard_Real aTol)
 150 {
 151   myTolTransition = aTol;
 152 }
 153 //=======================================================================
 154 //function : TolTransition
 155 //purpose  :
 156 //=======================================================================
 157   Standard_Real IntPatch_ALineToWLine::TolTransition()const
 158 {
 159   return myTolTransition;
 160 }
 161 //=======================================================================
 162 //function : SetTolParam
 163 //purpose  :
 164 //=======================================================================
 165   void IntPatch_ALineToWLine::SetTolParam(const Standard_Real aTolParam)
 166 {
 167   myTolParam = aTolParam;
 168 }
 169 //=======================================================================
 170 //function : TolParam
 171 //purpose  :
 172 //=======================================================================
 173   Standard_Real IntPatch_ALineToWLine::TolParam()const
 174 {
 175   return myTolParam;
 176 }
 177 //=======================================================================
 178 //function : SetTolOpenDomain
 179 //purpose  :
 180 //=======================================================================
 181   void IntPatch_ALineToWLine::SetTolOpenDomain(const Standard_Real aTol)
 182 {
 183   myTolOpenDomain = aTol;
 184 }
 185 //=======================================================================
 186 //function : TolOpenDomain
 187 //purpose  :
 188 //=======================================================================
 189   Standard_Real IntPatch_ALineToWLine::TolOpenDomain()const
 190 {
 191   return myTolOpenDomain;
 192 }
 193 //=======================================================================
 194 //function : SetConstantParameter
 195 //purpose  :
 196 //=======================================================================
 197   void IntPatch_ALineToWLine::SetConstantParameter() const
 198 {
 199 }
 200 //=======================================================================
 201 //function : SetUniformAbscissa
 202 //purpose  :
 203 //=======================================================================
 204   void IntPatch_ALineToWLine::SetUniformAbscissa() const
 205 {
 206 }
 207 //=======================================================================
 208 //function : SetUniformDeflection
 209 //purpose  :
 210 //=======================================================================
 211   void IntPatch_ALineToWLine::SetUniformDeflection() const
 212 {
 213 }
 214 //=======================================================================
 215 //function : MakeWLine
 216 //purpose  :
 217 //=======================================================================
 218   Handle(IntPatch_WLine) IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& aline) const
 219 {
 220   Standard_Boolean included;
 221   Standard_Real f = aline->FirstParameter(included);
 222   if(!included) {
 223     f+=myTolOpenDomain;
 224   }
 225   Standard_Real l = aline->LastParameter(included);
 226   if(!included) {
 227     l-=myTolOpenDomain;
 228   }
 229   return(MakeWLine(aline,f,l));
 230 }
 231
 232 //=======================================================================
 233 //function : MakeWLine
 234 //purpose  :
 235 //=======================================================================
 236   Handle(IntPatch_WLine) IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& aline,
 237                                                     const Standard_Real _firstparam,
 238                                                     const Standard_Real _lastparam) const
 239 {
 240   //
 241   Standard_Real dl, dumin, dumax, U, pv, pu1, pv1, pu2, pv2;
 242   Standard_Real firstparam, lastparam;
 243   Standard_Integer v, nbvtx;
 244   Standard_Boolean TriOk;
 245   //
 246   firstparam = _firstparam;
 247   lastparam  = _lastparam;
 248   // Pas Bon ...
 249   // nbpointsmax It is the field. ( =200. by default)
 250   dl = (lastparam - firstparam)/(Standard_Real)(nbpointsmax-1);
 251   dumin = 0.1*dl;
 252   dumax = 10 * dl;
 253   //
 254   nbvtx = aline->NbVertex();
 255   //
 256   TColStd_Array1OfReal paramvertex(1,Max(nbvtx,1)), newparamvertex(1,Max(nbvtx,1));
 257   //
 258   for(v = 1; v<=nbvtx; v++) {
 259     const IntPatch_Point& aVtx = aline->Vertex(v);
 260     pv=aVtx.ParameterOnLine();
 261     paramvertex(v)=pv;
 262     newparamvertex(v)=-1.;
 263   }
 264   //
 265   //-- Tri et Confusion des vertex proches
 266   do {
 267     TriOk = Standard_True;
 268     for(v=2; v<=nbvtx;v++) {
 269       if(paramvertex(v) < paramvertex(v-1)) {
 270         pv=paramvertex(v);
 271         paramvertex(v-1)=paramvertex(v);
 272         paramvertex(v)=pv;
 273         TriOk = Standard_False;
 274       }
 275     }
 276   }
 277   while(!TriOk);
 278   //
 279   for(v=2; v<=nbvtx;v++) {
 280     pv=paramvertex(v);
 281     pv1=paramvertex(v-1);
 282     if(pv-pv1 < myTolParam) {
 283       paramvertex(v)=pv1;
 284     }
 285   }
 286   //
 287   Standard_Integer t, nbpwline;
 288   Standard_Real u1,v1,u2,v2, anu1, anv1, anu2, anv2;
 289   Standard_Real aL, dl_sur_2;
 290   gp_Pnt Pnt3d, aPnt3d1;
 291   IntSurf_PntOn2S POn2S;
 292   Handle(IntSurf_LineOn2S) LinOn2S;
 293   //
 294   LinOn2S = new IntSurf_LineOn2S;
 295
 296   //// Modified by jgv, 17.09.09 for OCC21255 ////
 297   Standard_Real refpar = RealLast(), ref_u1 = 0., ref_u2 = 0.;
 298   if (nbvtx)
 299     {
 300       const IntPatch_Point& FirstVertex = aline->Vertex(1);
 301       refpar = FirstVertex.ParameterOnLine();
 302       FirstVertex.Parameters(ref_u1, v1, ref_u2, v2);
 303     }
 304   ////////////////////////////////////////////////
 305
 306   //-----------------------------------------------------
 307   //-- Estimation Grossiere de la longueur de la courbe
 308   //--
 309   aL=0.;
 310   dl_sur_2=0.5*dl;
 311   Pnt3d = aline->Value(firstparam);
 312   for(t=0, U=firstparam+dumin; U<lastparam; t++,U+=dl_sur_2) {
 313     aPnt3d1 = aline->Value(U);
 314     aL+=Pnt3d.Distance(aPnt3d1);
 315     Pnt3d=aPnt3d1;
 316   }
 317   //------------------------------------------------------
 318   aL/=t;
 319   aL*=2;
 320   //------------------------------------------------------
 321   //-- Calcul du premier point
 322   //--
 323   //------------------------------------------------------
 324   //-- On reajuste firstparam et lastparam au cas ou ces
 325   //-- valeurs sont tres proches des parametres de vertex
 326   if(nbvtx) {
 327     Standard_Real pvtxmin, pvtxmax;
 328     //
 329     pvtxmin = aline->Vertex(1).ParameterOnLine();
 330     pvtxmax = pvtxmin;
 331     //
 332     for(v=2; v<=nbvtx; v++) {
 333       pv=aline->Vertex(v).ParameterOnLine();
 334       if(pvtxmin>pv){
 335         pvtxmin = pv;
 336       }
 337       if(pvtxmax<pv){
 338         pvtxmax =pv;
 339       }
 340     }
 341     if(Abs(pvtxmin-firstparam)<myTolOpenDomain) {
 342       firstparam=pvtxmin;
 343     }
 344     if(Abs(pvtxmax-lastparam)<myTolOpenDomain) {
 345       lastparam=pvtxmax;
 346     }
 347   }
 348   //
 349   // First Point
 350   Pnt3d = aline->Value(firstparam);
 351   quad1.Parameters(Pnt3d,u1,v1);
 352   //
 353   RefineParameters(aline, firstparam, lastparam, firstparam, 1, quad1, 10.*myTol3D, u1,v1);
 354   //
 355   //
 356   quad2.Parameters(Pnt3d,u2,v2);
 357   //
 358   RefineParameters(aline, firstparam, lastparam, firstparam, 1, quad2, 10.*myTol3D, u2,v2);
 359   //
 360
 361   POn2S.SetValue(Pnt3d,u1,v1,u2,v2);
 362   anu1 = u1;
 363   anu2 = u2;
 364   LinOn2S->Add(POn2S);
 365   nbpwline = 1;
 366   U = firstparam;
 367   //-------------------------------------------------------
 368   //-- On detecte le cas : Point de debut == vertex
 369   //-- On affecte un parametre bidon qui sera recadre
 370   //-- dans l insertion des vertex dans la ligne
 371   //--
 372   for(v=1;v<=nbvtx;v++) {
 373     if(newparamvertex(v)<0.) {
 374       pv=paramvertex(v);
 375       if((pv>=U-2.0*myTolOpenDomain) && (pv<=U+2.0*myTolOpenDomain)) {
 376         if(pv-U > myTolParam) {
 377           newparamvertex(v) = 1.000001;
 378         }
 379         else if(U-pv>myTolParam) {
 380           newparamvertex(v) = 0.999999;
 381         }
 382       }
 383     }
 384   }
 385   //
 386   Standard_Real DeltaU;
 387   //
 388   // the loop
 389   DeltaU=0.;
 390   //// Modified by jgv, 17.09.09 for OCC21255 ////
 391   Standard_Boolean Corrected = Standard_False;
 392   ////////////////////////////////////////////////
 393   while(U<lastparam) {
 394     Standard_Integer NbCalculD1;
 395     Standard_Real UPourCalculD1, pvavant, pvapres;
 396     gp_Vec VPourCalculD1;
 397     gp_Pnt PPourCalculD1;
 398     //
 399     DeltaU=0.;
 400     NbCalculD1=0;
 401     UPourCalculD1=U;
 402     do {
 403       if(aline->D1(UPourCalculD1,PPourCalculD1,VPourCalculD1)) {
 404         Standard_Real NormV = VPourCalculD1.Magnitude();
 405         if(NormV > 1.0e-16) {
 406           DeltaU = aL/NormV;
 407         }
 408       }
 409       if(DeltaU < dumin)
 410         DeltaU = dumin;
 411       else if (DeltaU > dumax)
 412         DeltaU = dumax;
 413       UPourCalculD1+=dumin;
 414     }
 415     while((++NbCalculD1<10)&&(DeltaU==0.));
 416     //
 417     //OCC541(apo)->
 418     // static  //xft
 419     Standard_Real CurvDef = deflectionmax, AngDef = CurvDef;
 420     if(U+DeltaU < lastparam) {
 421       DefineDU(aline,U,DeltaU,CurvDef,AngDef);
 422     }
 423     //<-OCC451(apo)
 424     //
 425     if(DeltaU==0.){
 426       DeltaU = (dumin+dumax)*0.5;
 427     }
 428     //--------------------------------------------------------
 429     //-- On cherche a placer un minimum de ??? Points entre
 430     //-- 2 vertex
 431     //--------------------------------------------------------
 432     pvavant = firstparam;
 433     pvapres = lastparam;
 434     for(v=1;v<=nbvtx;v++) {
 435       pv=paramvertex(v);
 436       if(pv-U > myTolParam) {
 437         if(pvapres>pv) {
 438           pvapres = pv;
 439         }
 440       }
 441       else {
 442         if(U-pv > myTolParam) {
 443           if(pvavant<pv) {
 444             pvavant = pv;
 445           }
 446         }
 447       }
 448     }
 449     pvapres-=pvavant;
 450     if(pvapres < (10.*DeltaU)) {
 451       if(pvapres > (10.*dumin)) {
 452         DeltaU = pvapres*0.1;
 453       }
 454       else {
 455         DeltaU = dumin;
 456       }
 457     }
 458     //xf
 459     if (nbpwline==1 && nbvtx) {
 460       Standard_Real aUnext;
 461       //
 462       aUnext=U+DeltaU;
 463       pv=paramvertex(1);
 464       if (aUnext>pv){
 465         DeltaU=0.5*(pv-firstparam);
 466       }
 467     }
 468     //xt
 469     //--------------------------------------------------------
 470     //-- Calcul des nouveaux parametres sur les vertex
 471     //--
 472     for(v=1;v<=nbvtx;v++) {
 473       if(newparamvertex(v)<0.) {
 474         pv=paramvertex(v);
 475         if(pv>=U  && pv<(U+DeltaU) && (U+DeltaU<lastparam) ) {
 476           newparamvertex(v) = (Standard_Real)nbpwline + (pv-U)/DeltaU;
 477         }
 478       }
 479     }
 480
 481     //// Modified by jgv, 17.09.09 for OCC21255 ////
 482     if (!Corrected && U >= refpar)
 483       {
 484         CorrectFirstPartOfLine(LinOn2S, ref_u1, ref_u2, anu1, anu2);
 485         Corrected = Standard_True;
 486       }
 487     ////////////////////////////////////////////////
 488     U+=DeltaU;
 489     if(U < lastparam) {
 490       nbpwline++;
 491       Pnt3d = aline->Value(U);
 492       quad1.Parameters(Pnt3d,u1,v1);
 493       quad2.Parameters(Pnt3d,u2,v2);
 494       RecadreMemePeriode(u1,u2,anu1,anu2);
 495       anu1 = u1;
 496       anu2 = u2;
 497       POn2S.SetValue(Pnt3d,u1,v1,u2,v2);
 498       LinOn2S->Add(POn2S);
 499     }
 500   }//while(U<lastparam)
 501
 502   U-=DeltaU;
 503   for(v=1;v<=nbvtx;v++) {
 504     if(newparamvertex(v)<0.) {
 505       pv=paramvertex(v);
 506       if(pv <= lastparam+myTolOpenDomain) {
 507         if(lastparam-U) {
 508           newparamvertex(v) = (Standard_Real)nbpwline+(pv-U)/(lastparam-U);
 509         }
 510         else {
 511           newparamvertex(v) = nbpwline+1;
 512         }
 513       }
 514     }
 515   }
 516   //
 517   Pnt3d = aline->Value(lastparam);
 518   quad1.Parameters(Pnt3d,u1,v1);
 519   //
 520   RefineParameters(aline, firstparam, lastparam, lastparam, -1, quad1, 10.*myTol3D, u1,v1);
 521   //
 522   quad2.Parameters(Pnt3d,u2,v2);
 523   //
 524   RefineParameters(aline, firstparam, lastparam, lastparam, -1, quad2, 10.*myTol3D, u2,v2);
 525   //
 526   RecadreMemePeriode(u1,u2,anu1,anu2);
 527   POn2S.SetValue(Pnt3d,u1,v1,u2,v2);
 528   LinOn2S->Add(POn2S);
 529   nbpwline++;
 530
 531   //// Modified by jgv, 17.09.09 for OCC21255 ////
 532   if (!Corrected &&
 533       (lastparam >= refpar || refpar-lastparam < Precision::Confusion()))
 534     CorrectFirstPartOfLine(LinOn2S, ref_u1, ref_u2, anu1, anu2);
 535   ////////////////////////////////////////////////
 536
 537   //
 538   //-----------------------------------------------------------------
 539   //--  Calcul de la transition de la ligne sur les surfaces      ---
 540   //-----------------------------------------------------------------
 541   IntSurf_TypeTrans trans1,trans2;
 542   Standard_Integer indice1;
 543   Standard_Real dotcross;
 544   gp_Pnt aPP0, aPP1;
 545   //
 546   trans1=IntSurf_Undecided;
 547   trans2=IntSurf_Undecided;
 548   //
 549   indice1 = nbpwline/3;
 550   if(indice1<=2) {
 551     indice1 = 2;
 552   }
 553   //
 554   aPP1=LinOn2S->Value(indice1).Value();
 555   aPP0=LinOn2S->Value(indice1-1).Value();
 556   //
 557   gp_Vec tgvalid(aPP0, aPP1);
 558   gp_Vec aNQ1=quad1.Normale(aPP0);
 559   gp_Vec aNQ2=quad2.Normale(aPP0);
 560   //
 561   dotcross = tgvalid.DotCross(aNQ2, aNQ1);
 562   if (dotcross > myTolTransition) {
 563     trans1 = IntSurf_Out;
 564     trans2 = IntSurf_In;
 565   }
 566   else if(dotcross < -myTolTransition) {
 567     trans1 = IntSurf_In;
 568     trans2 = IntSurf_Out;
 569   }
 570   //-----------------------------------------------------------------
 571   //--              C r e a t i o n  d  e   la   W L i n e        ---
 572   //-----------------------------------------------------------------
 573   Handle(IntPatch_WLine) wline;
 574   //
 575   if(aline->TransitionOnS1() == IntSurf_Touch)  {
 576     Handle(IntPatch_WLine) wlinetemp = new IntPatch_WLine(LinOn2S,
 577                                               aline->IsTangent(),
 578                                               aline->SituationS1(),
 579                                               aline->SituationS2());
 580     wline = wlinetemp;
 581   }
 582   if(aline->TransitionOnS1() == IntSurf_Undecided)  {
 583     Handle(IntPatch_WLine) wlinetemp = new IntPatch_WLine(LinOn2S,
 584                                               aline->IsTangent());
 585     wline = wlinetemp;
 586   }
 587   else {
 588     Handle(IntPatch_WLine) wlinetemp = new IntPatch_WLine(LinOn2S,
 589                                               aline->IsTangent(),
 590                                               trans1, // aline->TransitionOnS1(),
 591                                               trans2);  //aline->TransitionOnS2());
 592     wline = wlinetemp;
 593   }
 594
 595   //-----------------------------------------------------------------
 596   //--         I n s e r t i o n   d  e s  v e r t e x            ---
 597   //-----------------------------------------------------------------
 598   TColStd_Array1OfInteger Redir(1,Max(nbvtx,1));
 599   //
 600   for(v=1;v<=nbvtx;v++) {
 601     Redir(v) = v;
 602   }
 603   //
 604   TriOk=Standard_True;
 605   Standard_Integer tamp;
 606   do {
 607     TriOk = Standard_True;
 608     for(v=2; v<=nbvtx;v++) {
 609       if(newparamvertex(Redir(v))<newparamvertex(Redir(v-1))) {
 610         tamp = Redir(v-1);
 611         Redir(v-1) = Redir(v);
 612         Redir(v) = tamp;
 613         TriOk = Standard_False;
 614       }
 615     }
 616   }
 617   while(!TriOk);
 618
 619   //-----------------------------------------------------------------
 620   //-- On detecte le cas ou un Vtx de param 1 ou nbpwline OnArc est double
 621   //-- Dans ce cas on supprime le vertex qui reference un arc deja reference
 622   //-- par un autre vertex
 623   //nbvtx = aline->NbVertex();
 624   Standard_Boolean APointHasBeenRemoved;
 625   //
 626   do {
 627     Standard_Boolean RemoveVtxo, RemoveVtx;
 628     Standard_Integer vo, voo;
 629     Standard_Real ponl, ponlo, aDist13, aDist23;
 630     //
 631     APointHasBeenRemoved = Standard_False;
 632     RemoveVtxo = Standard_False;
 633     RemoveVtx  = Standard_False;
 634     //
 635     for(v=1; v<=nbvtx   && !APointHasBeenRemoved; v++) {
 636       //
 637       if(newparamvertex(v)>=0.) {
 638         const IntPatch_Point& Vtx = aline->Vertex(v);
 639         ponl = Vtx.ParameterOnLine();
 640         const gp_Pnt& aP=Vtx.Value();
 641         //
 642         for(vo=1; vo<=nbvtx && !APointHasBeenRemoved; vo++) {
 643           if(v!=vo) {
 644             if(newparamvertex(vo)>=0.) {
 645               const IntPatch_Point& Vtxo = aline->Vertex(vo);
 646               ponlo = Vtxo.ParameterOnLine();
 647               const gp_Pnt& aPo=Vtxo.Value();
 648               //
 649               if(ponl-ponlo<myTolParam && ponlo-ponl<myTolParam) {
 650                 //
 651                 for(voo=1; voo<=nbvtx && !APointHasBeenRemoved; voo++) {
 652                   if(voo!=v && voo!=vo) {
 653                     if(newparamvertex(voo)>=0.) {
 654                       const IntPatch_Point& Vtxoo = aline->Vertex(voo);
 655                       const gp_Pnt& aPoo=Vtxoo.Value();
 656                       //
 657                       aDist13=aP.Distance(aPoo);
 658                       //
 659                       if(aDist13<=myTol3D) {
 660                         //-- 2 vertex de meme param  + un confondu geometriquement
 661                         if((Vtx.IsOnDomS1() == Vtxoo.IsOnDomS1()) &&
 662                            (Vtx.IsOnDomS2() == Vtxoo.IsOnDomS2()) ) {
 663
 664                           if(Vtx.IsOnDomS1()) {
 665                             if(Vtx.ParameterOnArc1()==Vtxoo.ParameterOnArc1()) {
 666                               if(SameCurve(Vtxoo.ArcOnS1(),Vtx.ArcOnS1())) { //-- param on s1 ?
 667                                 RemoveVtx = Standard_True;
 668                               }
 669                             }
 670                           }
 671                           if(Vtx.IsOnDomS2()) {
 672                             if(Vtx.ParameterOnArc2()==Vtxoo.ParameterOnArc2()) {
 673                               if(SameCurve(Vtxoo.ArcOnS2(),Vtx.ArcOnS2())) {
 674                                 RemoveVtx = Standard_True;
 675                               }
 676                             }
 677                             else {
 678                               RemoveVtx = Standard_False;
 679                             }
 680                           }
 681                         }
 682                         //
 683                         if(RemoveVtx==Standard_False) {
 684                           //
 685                           aDist23=aPo.Distance(aPoo);
 686                           //
 687                           if(aDist23<=myTol3D) {
 688                             //-- 2 vertex de meme param  + un confondu geometriquement
 689                             if((Vtxo.IsOnDomS1() == Vtxoo.IsOnDomS1())&&
 690                                (Vtxo.IsOnDomS2() == Vtxoo.IsOnDomS2()) ) {
 691                               if(Vtxo.IsOnDomS1()) {
 692                                 if(Vtxo.ParameterOnArc1()==Vtxoo.ParameterOnArc1()) {
 693                                   if(SameCurve(Vtxoo.ArcOnS1(),Vtxo.ArcOnS1())) { //-- param on s1 ?
 694                                     RemoveVtxo = Standard_True;
 695                                   }
 696                                 }
 697                               }
 698                               if(Vtxo.IsOnDomS2()) {
 699                                 if(Vtxo.ParameterOnArc2()==Vtxoo.ParameterOnArc2()) {
 700                                   if(SameCurve(Vtxoo.ArcOnS2(),Vtxo.ArcOnS2())) {
 701                                     RemoveVtxo = Standard_True;
 702                                   }
 703                                 }
 704                                 else {
 705                                   RemoveVtxo = Standard_False;
 706                                 }
 707                               }
 708                             }
 709                           }
 710                         } //--
 711                         if(RemoveVtx) {
 712                           newparamvertex(v) = -1.;
 713                           APointHasBeenRemoved = Standard_True;
 714                         }
 715                         else if(RemoveVtxo) {
 716                           newparamvertex(vo) = -1.;
 717                           APointHasBeenRemoved = Standard_True;
 718                         }
 719                       } //-- Pnt = Pntoo
 720                     } //-- voo!=v && voo!=vo
 721                   } //-- Fin boucle sur voo
 722                 }
 723               }
 724             }
 725           }
 726         }
 727       }
 728     }//for(v=1; v<=nbvtx   && !APointHasBeenRemoved; v++)
 729   }
 730   while(APointHasBeenRemoved);
 731   //
 732   Standard_Integer ParamVtxPrecedent, refpointonwline, aIndV;
 733   Standard_Real pvtx, approxparamonwline, aDst;
 734   Standard_Boolean bIsApex1, bIsApex2;
 735   //
 736   ParamVtxPrecedent = 0;
 737   //
 738   for(v=1; v<=nbvtx; v++) {
 739     aIndV=Redir(v);
 740     pv=paramvertex(v);// parameter from ALine
 741     pvtx = newparamvertex(aIndV);
 742     if(pvtx>=0. && (pvtx <= nbpwline+1)) {
 743       approxparamonwline=newparamvertex(aIndV);
 744       refpointonwline=1;
 745       //
 746       IntPatch_Point NewPoint = aline->Vertex(aIndV);
 747       //
 748       Pnt3d = NewPoint.Value();
 749       quad1.Parameters(Pnt3d, u1, v1);
 750       quad2.Parameters(Pnt3d, u2, v2);
 751       //-------------------------------------------------------
 752       //-- On recadre les parametres des vertex dans la bonne -
 753       //-- periode en recadrant avec le point le plus proche  -
 754       //-------------------------------------------------------
 755       if(approxparamonwline > nbpwline) {
 756         refpointonwline = nbpwline-1;
 757       }
 758       else if(approxparamonwline < 1) {
 759         refpointonwline = 1;
 760       }
 761       else {
 762         refpointonwline = (Standard_Integer)approxparamonwline;
 763       }
 764       //
 765       //
 766       const IntSurf_PntOn2S& aP2Sx=LinOn2S->Value(refpointonwline);
 767       aP2Sx.ParametersOnS1(anu1, anv1);
 768       aP2Sx.ParametersOnS2(anu2, anv2);
 769       //
 770       bIsApex1=IsApex(quad1, v1, myTol3D);
 771       bIsApex2=IsApex(quad2, v2, myTol3D);
 772       //
 773       //if (refpointonwline==1 || refpointonwline==nbpwline) {
 774       if (bIsApex1 || bIsApex2) {
 775         if (fabs(pv-firstparam)<myTolParam || fabs(pv-lastparam)<myTolParam) {
 776           // aline starts(ends) on vertex
 777           const gp_Pnt& aP1x=aP2Sx.Value();
 778           //
 779           aDst=aP1x.Distance(Pnt3d);
 780           if (aDst<10.*myTol3D) {
 781             u1=anu1;
 782             v1=anv1;
 783             u2=anu2;
 784             v2=anv2;
 785           }
 786         }
 787         //
 788         else {
 789           // aline goes through vertex
 790           if (bIsApex1) {
 791             u1=0.;
 792           }
 793           if (bIsApex2) {
 794             u2=0.;
 795           }
 796         }
 797       }
 798       //
 799       //
 800       if(v==1) {
 801         ParamVtxPrecedent=refpointonwline;
 802         RecadreMemePeriode(u1,u2,anu1,anu2);
 803         NewPoint.SetParameter(refpointonwline);
 804         //
 805         NewPoint.SetParameters(u1,v1,u2,v2);
 806         wline->AddVertex(NewPoint);
 807       }
 808       //
 809       else {
 810         if(ParamVtxPrecedent==refpointonwline) {
 811           //-- 2 vertex renseignent le meme point de la LineOn2S
 812           //-- On insere un nv point  =  vtx
 813           //-- On decale tous les vtx apres de 1
 814           RecadreMemePeriode(u1,u2,anu1,anu2);
 815           POn2S.SetValue(Pnt3d,u1,v1,u2,v2);
 816           LinOn2S->InsertBefore(refpointonwline+1, POn2S);
 817           nbpwline++;
 818           NewPoint.SetParameter(refpointonwline+1);
 819           NewPoint.SetParameters(u1,v1,u2,v2);
 820           wline->AddVertex(NewPoint);
 821           ParamVtxPrecedent = refpointonwline+1;
 822           //
 823           Standard_Integer vv=v+1;
 824           for(; vv<=nbvtx; vv++) {
 825             if(newparamvertex(Redir(vv))!=-1.) {
 826               newparamvertex(Redir(vv))=newparamvertex(Redir(vv))+1.;
 827             }
 828           }
 829         }
 830         //
 831         else {
 832           RecadreMemePeriode(u1,u2, anu1, anu2);
 833           NewPoint.SetParameter(refpointonwline);
 834           //
 835           NewPoint.SetParameters(u1, v1, u2, v2);
 836           wline->AddVertex(NewPoint);
 837           ParamVtxPrecedent = refpointonwline;
 838         }
 839       }
 840     }
 841   }
 842   //
 843   pu1=0.;
 844   pv1=0.;
 845   pu2=0.;
 846   pv2=0.;
 847   //
 848   switch(quad1.TypeQuadric()) {
 849   case GeomAbs_Cylinder:
 850   case GeomAbs_Cone:
 851   case GeomAbs_Sphere:
 852     pu1=M_PI+M_PI;
 853     break;
 854   case GeomAbs_Torus:
 855     pu1=pv1=M_PI+M_PI;
 856     break;
 857   default:
 858     break;
 859   }
 860   switch(quad2.TypeQuadric()) {
 861   case GeomAbs_Cylinder:
 862   case GeomAbs_Cone:
 863   case GeomAbs_Sphere:
 864     pu2=M_PI+M_PI;
 865     break;
 866   case GeomAbs_Torus:
 867     pu2=pv2=M_PI+M_PI;
 868     break;
 869   default:
 870     break;
 871   }
 872   wline->SetPeriod(pu1,pv1,pu2,pv2);
 873
 874   wline->ComputeVertexParameters(myTol3D);
 875   return(wline);
 876 }
 877 //=======================================================================
 878 //function : DefineDU
 879 //purpose  :
 880 //=======================================================================
 881 gp_Pnt DefineDU(const Handle(IntPatch_ALine)& aline,
 882                 const Standard_Real U,
 883                 Standard_Real& DU,
 884                 const Standard_Real CurvDef,
 885                 const Standard_Real AngDef)
 886 {
 887   gp_Pnt P1 = aline->Value(U), P2, P3;
 888   gp_Vec V13, V12, V23;
 889   Standard_Real dU = DU/2.0, curvDef, angDef, m1, m2, m3;
 890   for(;;) { //According to class TangentialDeflection from GCPnts
 891     P2=aline->Value(U+dU);  P3=aline->Value(U+DU);
 892     V13 = P3.XYZ().Subtracted(P1.XYZ());  m1 = V13.Magnitude();
 893     V12 = P2.XYZ().Subtracted(P1.XYZ());  m2 = V12.Magnitude();
 894     V23 = P3.XYZ().Subtracted(P2.XYZ());  m3 = V23.Magnitude();
 895     if(m1 < CurvDef || m2 < CurvDef || m3 < CurvDef) break;
 896     curvDef = Abs(V13.Angle(V12));
 897     angDef = Abs(V13.Angle(V23));
 898     if(curvDef < CurvDef && angDef < AngDef) break;
 899     DU = dU;  dU /= 2.0;
 900   }
 901   return P3;
 902 }
 903 //=======================================================================
 904 //function : SameCurve
 905 //purpose  :
 906 //=======================================================================
 907 Standard_Boolean SameCurve(const Handle_Adaptor2d_HCurve2d& C1,const Handle_Adaptor2d_HCurve2d& C2)
 908 {
 909   Standard_Real C1f = C1->FirstParameter();
 910   Standard_Real C2f = C2->FirstParameter();
 911   if(C1f!=C2f) return(Standard_False);
 912   Standard_Real C1l = C1->LastParameter();
 913   Standard_Real C2l = C2->LastParameter();
 914   if(C1l!=C2l) return(Standard_False);
 915   Standard_Real u=0.3*C1f+0.7*C1l;
 916   gp_Pnt2d P1 = C1->Value(u);
 917   gp_Pnt2d P2 = C2->Value(u);
 918   if(P1.X()!=P2.X()) return(Standard_False);
 919   if(P1.Y()!=P2.Y()) return(Standard_False);
 920   return(Standard_True);
 921 }
 922 //=======================================================================
 923 //function : RecadreMemePeriode
 924 //purpose  :
 925 //=======================================================================
 926 void RecadreMemePeriode(Standard_Real& u1,
 927                         Standard_Real& u2,
 928                         const Standard_Real anu1,
 929                         const Standard_Real anu2)
 930 {
 931   //
 932   while(anu1-u1 > 5.0) {
 933     u1+=M_PI+M_PI;
 934   }
 935   while(u1-anu1 > 5.0) {
 936     //
 937     /*
 938     if (!bBothCylinders) {//cfe900/H6
 939       // this check on Cylinder/Cylinder intersection is probably
 940       // because of pbs with ALine on it.
 941       // For other Quadrics there was no pbs found during
 942       // grid tests.
 943       // (it is necessary to see ...IntCyCy(...) for details).
 944       //
 945       // In any case the pb does not deal with apex problem.
 946       //
 947       if (u1-M_PI-M_PI<0.) {
 948         break;
 949       }
 950     }
 951     */
 952     //
 953     u1-=M_PI+M_PI;
 954   }
 955   while(anu2-u2 > 5.0) {
 956     u2+=M_PI+M_PI;
 957   }
 958   while(u2-anu2 > 5.0) {
 959     //
 960     /*
 961     if (!bBothCylinders) {//cfe900/H6
 962       if (u2-M_PI-M_PI<0.) {
 963         break;
 964       }
 965     }
 966     */
 967     //
 968     u2-=M_PI+M_PI;
 969   }
 970 }
 971
 972 //=======================================================================
 973 //function : CorrectFirstPartOfLine
 974 //purpose  :
 975 //=======================================================================
 976 void CorrectFirstPartOfLine(Handle(IntSurf_LineOn2S)& LinOn2S,
 977                             const Standard_Real ref_u1,
 978                             const Standard_Real ref_u2,
 979                             Standard_Real& new_u1,
 980                             Standard_Real& new_u2)
 981 {
 982   Standard_Integer nbp = LinOn2S->NbPoints();
 983   Standard_Real u1, v1, u2, v2, OffsetOnS1, OffsetOnS2;
 984
 985   IntSurf_PntOn2S aPoint = LinOn2S->Value(nbp);
 986   aPoint.Parameters(u1, v1, u2, v2);
 987
 988   new_u1 = u1;
 989   new_u2 = u2;
 990   RecadreMemePeriode(new_u1, new_u2, ref_u1, ref_u2);
 991   OffsetOnS1 = new_u1 - u1;
 992   OffsetOnS2 = new_u2 - u2;
 993   if (Abs(OffsetOnS1) > 1. || Abs(OffsetOnS2) > 1.) //recadre on n*2*PI is done
 994     {
 995       Standard_Integer i;
 996       for (i = 1; i <= nbp; i++)
 997         {
 998           aPoint = LinOn2S->Value(i);
 999           aPoint.Parameters(u1, v1, u2, v2);
1000           LinOn2S->SetUV( i, Standard_True,  u1 + OffsetOnS1, v1 );
1001           LinOn2S->SetUV( i, Standard_False, u2 + OffsetOnS2, v2 );
1002         }
1003     }
1004 }
1005
1006 //
1007 //=======================================================================
1008 //function : IsApex
1009 //purpose  :
1010 //=======================================================================
1011 Standard_Boolean IsApex(const IntSurf_Quadric& aQuadric,
1012                         const Standard_Real aVx,
1013                         const Standard_Real aTol3D)
1014
1015 {
1016   Standard_Boolean bFlag;
1017   Standard_Real aHalfPi, aEpsilon;
1018   GeomAbs_SurfaceType aType;
1019   //
1020   bFlag=Standard_False;
1021   //
1022   aType=aQuadric.TypeQuadric();
1023   if (!(aType==GeomAbs_Cone || aType==GeomAbs_Sphere)) {
1024     return bFlag;
1025   }
1026   //
1027   aEpsilon=Epsilon(10.);//1.77e-15
1028   //
1029   // apex on the Sphere
1030   if(aType==GeomAbs_Sphere) {
1031     aHalfPi=0.5*M_PI;
1032     if (fabs(aVx-aHalfPi)<aEpsilon) {
1033       bFlag=!bFlag;
1034     }
1035     else if (fabs(aVx+aHalfPi)<aEpsilon){
1036       bFlag=!bFlag;
1037     }
1038   }
1039   //
1040   // apex on the Cone
1041   else if(aType==GeomAbs_Cone) {
1042     Standard_Real aDst;
1043     gp_Pnt aPap, aPx;
1044     gp_Cone aCone;
1045     //
1046     aCone=aQuadric.Cone();
1047     aPap=aCone.Apex();
1048     aPx=aQuadric.Value(0.,aVx);
1049     //
1050     aDst=aPx.Distance(aPap);
1051     if(aDst<aTol3D) {
1052       bFlag=!bFlag;
1053     }
1054   }
1055   return bFlag;
1056 }
1057 //=======================================================================
1058 //function : RefineParameters
1059 //purpose  :
1060 //=======================================================================
1061 void RefineParameters(const Handle(IntPatch_ALine)& aALine,
1062                       const Standard_Real aTb,
1063                       const Standard_Real aTe,
1064                       const Standard_Real aTx,
1065                       const Standard_Integer iDir,
1066                       const IntSurf_Quadric& aQuadric,
1067                       const Standard_Real aTol3D,
1068                       Standard_Real& aUx,
1069                       Standard_Real& aVx)
1070 {
1071   GeomAbs_SurfaceType aType;
1072   //
1073   aType=aQuadric.TypeQuadric();
1074   if (!(aType==GeomAbs_Cone || aType==GeomAbs_Sphere)) {
1075     return;
1076   }
1077   //
1078   Standard_Boolean bIsDone, bIsEmpty, bParallel, bFound;
1079   Standard_Integer aNbPoints = 0;
1080   Standard_Real aHalfPi, aEpsilon, aLimV, dT, aT1, aT2 = 0., aEpsT;
1081   Standard_Real aU1, aV1, aU2, aV2;
1082   gp_Pnt aP1, aP2, aPx;
1083   gp_Pnt2d  aP2D1, aP2D2, aPLim(0., 0.);
1084   gp_Vec2d aVLim(1., 0.);
1085   gp_Lin2d aLLim;
1086   IntAna2d_AnaIntersection aAI;
1087   //
1088   aEpsilon=Epsilon(10.);//1.77e-15
1089   aEpsT=0.0001;
1090   aLLim.SetDirection(aVLim);
1091   //
1092   // apex on the Cone
1093   if(aType==GeomAbs_Cone) {
1094     Standard_Real aDst;
1095     gp_Pnt aPap;
1096     gp_Cone aCone;
1097     //
1098     aCone=aQuadric.Cone();
1099     aPap=aCone.Apex();
1100     //aPx=aQuadric.Value(0.,aVx);
1101     aPx=aALine->Value(aTx);
1102     //
1103     aDst=aPx.Distance(aPap);
1104     if(aDst>aTol3D) {// nothing to do
1105       return;
1106     }
1107     //
1108     aPLim.SetY(aVx);
1109     aLLim.SetLocation(aPLim);
1110     //
1111   }
1112   //
1113   // apex on the Sphere
1114   if(aType==GeomAbs_Sphere) {
1115     aHalfPi=0.5*M_PI;
1116     //
1117     if (fabs(aVx-aHalfPi)<aEpsilon) {
1118       aLimV=aHalfPi;
1119     }
1120     else if (fabs(aVx+aHalfPi)<aEpsilon){
1121       aLimV=-aHalfPi;
1122     }
1123     else {
1124       //Check: aUx must be 0 or 2*pi
1125       if(fabs(aUx) < aEpsilon || fabs(aUx - 2.*M_PI) < aEpsilon) {
1126         //aUx = 0 or 2*pi, but may be it must be 2*pi or 0?
1127         bFound=FindNearParameter(aALine, aTx, iDir, aTol3D, aT1);
1128         if(!bFound) {
1129           dT=aEpsT*(aTe-aTb);
1130           if (iDir<0) {
1131             dT=-dT;
1132           }
1133           aT1=aTx+dT;
1134         }
1135
1136         aP1=aALine->Value(aT1);
1137         aQuadric.Parameters(aP1, aU1, aV1);
1138
1139         if(fabs(aU1) > fabs(aU1 - 2.*M_PI)) {
1140           aUx = 2.*M_PI;
1141         }
1142         else {
1143           aUx = 0.;
1144         }
1145       }
1146
1147       return;
1148     }
1149     //
1150     aPLim.SetY(aLimV);
1151     aLLim.SetLocation(aPLim);
1152   }
1153   //
1154   // aT1, aT2
1155   //
1156   // Try to find aT1, aT2 taking into acc 3D Tolerance
1157   bFound=FindNearParameter(aALine, aTx, iDir, aTol3D, aT1);
1158   if (bFound) {
1159     bFound=FindNearParameter(aALine, aT1, iDir, aTol3D, aT2);
1160   }
1161   if (!bFound) {
1162     // Assign aT1, aT2 by some values
1163     dT=aEpsT*(aTe-aTb);
1164     if (iDir<0) {
1165       dT=-dT;
1166     }
1167     aT1=aTx+dT;
1168     aT2=aT1+dT;
1169   }
1170   //
1171   aP1=aALine->Value(aT1);
1172   aQuadric.Parameters(aP1, aU1, aV1);
1173   aP2D1.SetCoord(aU1, aV1);
1174   //
1175   aP2=aALine->Value(aT2);
1176   aQuadric.Parameters(aP2, aU2, aV2);
1177   aP2D2.SetCoord(aU2, aV2);
1178   //
1179   gp_Vec2d aV12(aP2D1, aP2D2);
1180
1181   if(aV12.SquareMagnitude() <= aEpsilon) {
1182     return;
1183   }
1184   gp_Lin2d aL12(aP2D1, aV12);
1185   //
1186   aAI.Perform(aLLim, aL12);
1187   bIsDone=aAI.IsDone();
1188   if (!bIsDone) {
1189     return;
1190   }
1191   bIsEmpty=aAI.IsEmpty();
1192   if (bIsEmpty) {
1193     return;
1194   }
1195   aNbPoints=aAI.NbPoints();
1196   if (!aNbPoints) {
1197     return;
1198   }
1199   bParallel=aAI.ParallelElements();
1200   if (bParallel) {
1201     return;
1202   }
1203   const IntAna2d_IntPoint& aIAPnt=aAI.Point(1);
1204   const gp_Pnt2d& aP2D=aIAPnt.Value();
1205   aUx=aP2D.X();
1206 }
1207 //=======================================================================
1208 //function : FindNearParameter
1209 //purpose  :
1210 //=======================================================================
1211 Standard_Boolean FindNearParameter(const Handle(IntPatch_ALine)& aALine,
1212                                    const Standard_Real aTx,
1213                                    const Standard_Integer iDir,
1214                                    const Standard_Real aTol3D,
1215                                    Standard_Real& aT1)
1216 {
1217   Standard_Boolean bFound;
1218   Standard_Real aX, aY, aZ;
1219   gp_Pnt aPx, aP1;
1220   gp_Vec aVx;
1221   //
1222   aT1=0.;
1223   //
1224   bFound=aALine->D1(aTx, aPx, aVx);
1225   if(!bFound){
1226     return bFound;
1227   }
1228   gp_Dir aDx(aVx);
1229   if (iDir<0) {
1230     aDx.Reverse();
1231   }
1232   aX=aPx.X()+aDx.X()*aTol3D;
1233   aY=aPx.Y()+aDx.Y()*aTol3D;
1234   aZ=aPx.Z()+aDx.Z()*aTol3D;
1235   aP1.SetCoord(aX, aY, aZ);
1236   //
1237   bFound=aALine->FindParameter(aP1, aT1);
1238   //
1239   return bFound;
1240 }
1241