src/GeometryTest/GeometryTest_CurveCommands.cxx

   1 // Created on: 1993-08-12
   2 // Created by: Bruno DUMORTIER
   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 under
   9 // the terms of the GNU Lesser General Public License 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 // 09/06/97 : JPI : suppression des commandes redondantes suite a la creation de GeomliteTest
  18
  19 #include <GeometryTest.hxx>
  20 #include <Draw_Appli.hxx>
  21 #include <DrawTrSurf.hxx>
  22 #include <DrawTrSurf_BSplineCurve.hxx>
  23 #include <DrawTrSurf_BSplineCurve2d.hxx>
  24 #include <Draw_Marker3D.hxx>
  25 #include <Draw_Marker2D.hxx>
  26 #include <Draw.hxx>
  27 #include <Draw_Interpretor.hxx>
  28 #include <Draw_Color.hxx>
  29
  30 #include <GeomAPI.hxx>
  31 #include <GeomAPI_IntCS.hxx>
  32 #include <GeomAPI_IntSS.hxx>
  33
  34 //#include <GeomLProp.hxx>
  35 #include <GeomProjLib.hxx>
  36 #include <BSplCLib.hxx>
  37
  38 #include <gp.hxx>
  39 #include <gp_Pln.hxx>
  40
  41 #include <Geom_Line.hxx>
  42 #include <Geom_Ellipse.hxx>
  43 #include <Geom2d_Line.hxx>
  44 #include <Geom2d_Circle.hxx>
  45 #include <Geom2d_Ellipse.hxx>
  46 #include <Geom2d_Parabola.hxx>
  47 #include <Geom2d_Hyperbola.hxx>
  48 #include <Geom2d_BSplineCurve.hxx>
  49 #include <Geom2d_Curve.hxx>
  50
  51 #include <GccAna_Lin2dBisec.hxx>
  52 #include <GccAna_Circ2dBisec.hxx>
  53 #include <GccAna_CircLin2dBisec.hxx>
  54 #include <GccAna_CircPnt2dBisec.hxx>
  55 #include <GccAna_LinPnt2dBisec.hxx>
  56 #include <GccAna_Pnt2dBisec.hxx>
  57 #include <GccInt_Bisec.hxx>
  58 #include <GccInt_IType.hxx>
  59
  60 #include <Geom_Plane.hxx>
  61 #include <Geom_Curve.hxx>
  62
  63 #include <Law_BSpline.hxx>
  64
  65 #include <TColgp_Array1OfPnt.hxx>
  66 #include <TColgp_Array1OfPnt2d.hxx>
  67 #include <TColStd_Array1OfReal.hxx>
  68 #include <TColStd_Array1OfInteger.hxx>
  69
  70 #include <Adaptor3d_Curve.hxx>
  71
  72 #include <GeomAdaptor_Surface.hxx>
  73
  74 #include <ProjLib_HCompProjectedCurve.hxx>
  75 #include <Precision.hxx>
  76 #include <Message.hxx>
  77
  78 #include <Geom_Surface.hxx>
  79 #include <stdio.h>
  80 #include <Geom_BSplineCurve.hxx>
  81 #include <GCPnts_QuasiUniformDeflection.hxx>
  82 #include <GCPnts_TangentialDeflection.hxx>
  83 #include <GCPnts_DistFunction.hxx>
  84 #include <gce_MakeLin.hxx>
  85 #include <TColStd_Array1OfBoolean.hxx>
  86 #include <Adaptor3d_TopolTool.hxx>
  87 #include <TColgp_Array2OfPnt.hxx>
  88 #include <Geom_BSplineSurface.hxx>
  89 #include <DrawTrSurf_BSplineSurface.hxx>
  90
  91 //epa test
  92 #include <BRepBuilderAPI_MakeEdge.hxx>
  93 #include <AIS_Shape.hxx>
  94 #include <TopoDS.hxx>
  95 #include <BRepAdaptor_CompCurve.hxx>
  96 #include <GeomLProp_CLProps.hxx>
  97 #include <GCPnts_AbscissaPoint.hxx>
  98 #include <GCPnts_UniformAbscissa.hxx>
  99 #include <DBRep.hxx>
 100
 101 #ifdef _WIN32
 102 Standard_IMPORT Draw_Viewer dout;
 103 #endif
 104
 105 //=======================================================================
 106 //function : polecurve2d
 107 //purpose  :
 108 //=======================================================================
 109
 110 static Standard_Integer polelaw (Draw_Interpretor& , Standard_Integer n, const char** a)
 111 {
 112   Standard_Integer k,
 113   jj,
 114   qq,
 115   i;
 116
 117
 118   if (n < 3) return 1;
 119   Standard_Boolean periodic = Standard_False ;
 120   Standard_Integer deg = Draw::Atoi(a[2]);
 121   Standard_Integer nbk = Draw::Atoi(a[3]);
 122
 123   TColStd_Array1OfReal    knots(1, nbk);
 124   TColStd_Array1OfInteger mults(1, nbk);
 125   k = 4;
 126   Standard_Integer Sigma = 0;
 127   for (i = 1; i<=nbk; i++) {
 128     knots( i) = Draw::Atof(a[k]);
 129     k++;
 130     mults( i) = Draw::Atoi(a[k]);
 131     Sigma += mults(i);
 132     k++;
 133   }
 134
 135   Standard_Integer np;
 136   np = Sigma - deg  -1;
 137   TColStd_Array1OfReal flat_knots(1, Sigma) ;
 138   jj = 1 ;
 139   for (i = 1 ; i <= nbk ; i++) {
 140     for(qq = 1 ; qq <= mults(i) ; qq++) {
 141       flat_knots(jj) = knots(i) ;
 142       jj ++ ;
 143     }
 144   }
 145
 146   TColgp_Array1OfPnt2d poles  (1, np);
 147   TColStd_Array1OfReal schoenberg_points(1,np) ;
 148   BSplCLib::BuildSchoenbergPoints(deg,
 149                                   flat_knots,
 150                                   schoenberg_points) ;
 151   for (i = 1; i <= np; i++) {
 152     poles(i).SetCoord(schoenberg_points(i),Draw::Atof(a[k]));
 153     k++;
 154   }
 155
 156   Handle(Geom2d_BSplineCurve) result =
 157     new Geom2d_BSplineCurve(poles, knots, mults, deg, periodic);
 158   DrawTrSurf::Set(a[1],result);
 159
 160
 161   return 0;
 162 }
 163 //=======================================================================
 164 //function : to2d
 165 //purpose  :
 166 //=======================================================================
 167
 168 static Standard_Integer to2d (Draw_Interpretor& , Standard_Integer n, const char** a)
 169 {
 170   if (n < 3) return 1;
 171
 172   // get the curve
 173   Handle(Geom_Curve) C = DrawTrSurf::GetCurve(a[2]);
 174   if (C.IsNull())
 175     return 1;
 176
 177   Handle(Geom_Surface) S;
 178   if (n >= 4) {
 179     S = DrawTrSurf::GetSurface(a[3]);
 180     if (S.IsNull()) return 1;
 181   }
 182   else
 183     S = new Geom_Plane(gp::XOY());
 184
 185   Handle(Geom_Plane) P = Handle(Geom_Plane)::DownCast(S);
 186   if (P.IsNull()) return 1;
 187   Handle(Geom2d_Curve) r = GeomAPI::To2d(C,P->Pln());
 188   DrawTrSurf::Set(a[1],r);
 189   return 0;
 190 }
 191
 192 //=======================================================================
 193 //function : to3d
 194 //purpose  :
 195 //=======================================================================
 196
 197 static Standard_Integer to3d (Draw_Interpretor& , Standard_Integer n, const char** a)
 198 {
 199   if (n < 3) return 1;
 200
 201   Handle(Geom2d_Curve) C = DrawTrSurf::GetCurve2d(a[2]);
 202   if (C.IsNull()) return 1;
 203
 204   Handle(Geom_Surface) S;
 205   if (n >= 4) {
 206     S = DrawTrSurf::GetSurface(a[3]);
 207     if (S.IsNull()) return 1;
 208   }
 209   else
 210     S = new Geom_Plane(gp::XOY());
 211
 212   Handle(Geom_Plane) P = Handle(Geom_Plane)::DownCast(S);
 213   if (P.IsNull()) return 1;
 214   Handle(Geom_Curve) r = GeomAPI::To3d(C,P->Pln());
 215
 216   DrawTrSurf::Set(a[1],r);
 217   return 0;
 218 }
 219
 220 //=======================================================================
 221 //function : gproject
 222 //purpose  :
 223 //=======================================================================
 224
 225
 226 static Standard_Integer gproject(Draw_Interpretor& di, Standard_Integer n, const char** a)
 227 {
 228   TCollection_AsciiString newname;
 229   TCollection_AsciiString newname1;
 230
 231   if (n < 4)
 232   {
 233     di << "gproject waits 3 or more arguments\n";
 234     return 1;
 235   }
 236
 237   TCollection_AsciiString name = a[1];
 238
 239   Handle(Geom_Curve) Cur = DrawTrSurf::GetCurve(a[2]);
 240   Handle(Geom_Surface) Sur = DrawTrSurf::GetSurface(a[3]);
 241   if (Cur.IsNull() || Sur.IsNull()) return 1;
 242
 243   Handle(GeomAdaptor_Curve) hcur = new GeomAdaptor_Curve(Cur);
 244   Handle(GeomAdaptor_Surface) hsur = new GeomAdaptor_Surface(Sur);
 245
 246   Standard_Integer index = 4;
 247   Standard_Real aTol3d = 1.e-6;
 248   Standard_Real aMaxDist = -1.0;
 249
 250   if (n > 4 && a[4][0] != '-')
 251   {
 252     aTol3d = Draw::Atof(a[4]);
 253     index = 5;
 254
 255     if (n > 5 && a[5][0] != '-')
 256     {
 257       aMaxDist = Draw::Atof(a[5]);
 258       index = 6;
 259     }
 260   }
 261
 262   Handle(ProjLib_HCompProjectedCurve) HProjector = new ProjLib_HCompProjectedCurve(aTol3d, hsur, hcur, aMaxDist);
 263   ProjLib_CompProjectedCurve& Projector = *HProjector;
 264
 265   GeomAbs_Shape aContinuity = GeomAbs_C2;
 266   Standard_Integer aMaxDegree, aMaxSeg;
 267   Standard_Boolean aProj2d;
 268   Standard_Boolean aProj3d;
 269
 270   while (index + 1 < n)
 271   {
 272     if (a[index][0] != '-') return 1;
 273
 274     if (a[index][1] == 'c')
 275     {
 276       Standard_CString aContinuityName = a[index + 1];
 277       if (!strcmp(aContinuityName, "C0"))
 278       {
 279         aContinuity = GeomAbs_C0;
 280       }
 281       else if (!strcmp(aContinuityName, "C1"))
 282       {
 283         aContinuity = GeomAbs_C1;
 284       }
 285       else if (!strcmp(aContinuityName, "C2"))
 286       {
 287         aContinuity = GeomAbs_C2;
 288       }
 289
 290       Projector.SetContinuity(aContinuity);
 291     }
 292     else if (a[index][1] == 'd')
 293     {
 294       aMaxDegree = Draw::Atoi(a[index + 1]);
 295       aMaxDegree = aMaxDegree > 25 ? 25 : aMaxDegree;
 296       Projector.SetMaxDegree(aMaxDegree);
 297     }
 298     else if (a[index][1] == 's')
 299     {
 300       aMaxSeg = Draw::Atoi(a[index + 1]);
 301       Projector.SetMaxSeg(aMaxSeg);
 302     }
 303     else if (!strcmp(a[index], "-2d"))
 304     {
 305       aProj2d = Draw::Atoi(a[index + 1]) > 0 ? Standard_True : Standard_False;
 306       Projector.SetProj2d(aProj2d);
 307     }
 308     else if (!strcmp(a[index], "-3d"))
 309     {
 310       aProj3d = Draw::Atoi(a[index + 1]) > 0 ? Standard_True : Standard_False;
 311       Projector.SetProj3d(aProj3d);
 312     }
 313
 314     index += 2;
 315   }
 316
 317   Projector.Perform();
 318
 319   for (Standard_Integer k = 1; k <= Projector.NbCurves(); k++) {
 320     newname  = name +   "_" + TCollection_AsciiString(k);
 321     newname1 = name + "2d_" + TCollection_AsciiString(k);
 322
 323     if (Projector.ResultIsPoint(k))
 324     {
 325       if (Projector.GetProj2d())
 326       {
 327         DrawTrSurf::Set(newname1.ToCString(), Projector.GetResult2dP(k));
 328         di << newname1 << " is pcurve\n";
 329       }
 330       if (Projector.GetProj3d())
 331       {
 332         DrawTrSurf::Set(newname.ToCString(), Projector.GetResult3dP(k));
 333         di << newname << " is 3d projected curve\n";
 334       }
 335     }
 336     else {
 337       if (Projector.GetProj2d())
 338       {
 339         DrawTrSurf::Set(newname1.ToCString(), Projector.GetResult2dC(k));
 340
 341         di << newname1 << " is pcurve\n";
 342         di << " Tolerance reached in 2d is " << Projector.GetResult2dUApproxError(k)
 343             << ";  " << Projector.GetResult2dVApproxError(k) << "\n";
 344       }
 345       if (Projector.GetProj3d())
 346       {
 347         DrawTrSurf::Set(newname.ToCString(), Projector.GetResult3dC(k));
 348
 349         di << newname << " is 3d projected curve\n";
 350         di << " Tolerance reached in 3d is " << Projector.GetResult3dApproxError(k) << "\n";
 351       }
 352     }
 353   }
 354   return 0;
 355 }
 356
 357 //=======================================================================
 358 //function : project
 359 //purpose  :
 360 //=======================================================================
 361
 362 static Standard_Integer project (Draw_Interpretor& di,
 363                                  Standard_Integer n, const char** a)
 364 {
 365   if ( n == 1) {
 366
 367     di << "project result2d c3d surf [-e p] [-v n] [-t tol]\n";
 368     di << "   -e p   : extent the surface of <p>%\n";
 369     di << "   -v n   : verify the projection at <n> points.\n";
 370     di << "   -t tol : set the tolerance for approximation\n";
 371     return 0;
 372   }
 373
 374   if (n < 4) return 1;
 375   Handle(Geom_Surface) GS = DrawTrSurf::GetSurface(a[3]);
 376   if (GS.IsNull()) return 1;
 377
 378   Handle(Geom_Curve) GC = DrawTrSurf::GetCurve(a[2]);
 379   if (GC.IsNull()) return 1;
 380
 381   Standard_Real tolerance = Precision::Confusion() ;
 382
 383   Standard_Real U1,U2,V1,V2;
 384   GS->Bounds(U1,U2,V1,V2);
 385
 386   Standard_Boolean Verif = Standard_False;
 387   Standard_Integer NbPoints=0;
 388
 389   Standard_Integer index = 4;
 390   while ( index+1 < n) {
 391     if ( a[index][0] != '-') return 1;
 392
 393     if ( a[index][1] == 'e') {
 394       Standard_Real p = Draw::Atof(a[index+1]);
 395       Standard_Real dU = p * (U2 - U1) / 100.;
 396       Standard_Real dV = p * (V2 - V1) / 100.;
 397       U1 -= dU; U2 += dU; V1 -= dV; V2 += dV;
 398     }
 399     else if ( a[index][1] == 'v') {
 400       Verif = Standard_True;
 401       NbPoints = Draw::Atoi(a[index+1]);
 402     }
 403     else if ( a[index][1] == 't') {
 404       tolerance = Draw::Atof(a[index+1]);
 405     }
 406     index += 2;
 407   }
 408
 409   Handle(Geom2d_Curve) G2d =
 410     GeomProjLib::Curve2d(GC, GS, U1, U2, V1, V2, tolerance);
 411
 412   if ( G2d.IsNull() ) {
 413     di << "\nProjection Failed\n";
 414     return 1;
 415   }
 416   else {
 417     DrawTrSurf::Set(a[1],G2d);
 418   }
 419   if ( Verif) {  // verify the projection on n points
 420     if ( NbPoints <= 0) {
 421       di << " n must be positive\n";
 422       return 0;
 423     }
 424     gp_Pnt P1,P2;
 425     gp_Pnt2d P2d;
 426
 427     Standard_Real U, dU;
 428     Standard_Real Dist,DistMax = -1.;
 429     U1 = GC->FirstParameter();
 430     U2 = GC->LastParameter();
 431     dU = ( U2 - U1) / (NbPoints + 1);
 432     for ( Standard_Integer i = 0 ; i <= NbPoints +1; i++) {
 433       U = U1 + i *dU;
 434       P1 = GC->Value(U);
 435       P2d = G2d->Value(U);
 436       P2 = GS->Value(P2d.X(), P2d.Y());
 437       Dist = P1.Distance(P2);
 438       di << " Parameter = " << U << "\tDistance = " << Dist << "\n";
 439       if ( Dist > DistMax) DistMax = Dist;
 440     }
 441     di << " **** Distance Maximale : " << DistMax << "\n";
 442   }
 443
 444   return 0;
 445 }
 446
 447 //=======================================================================
 448 //function : projonplane
 449 //purpose  :
 450 //=======================================================================
 451
 452 Standard_Integer projonplane(Draw_Interpretor& di,
 453                              Standard_Integer n, const char** a)
 454 {
 455   if ( n < 4 ) return 1;
 456
 457   Handle(Geom_Surface) S = DrawTrSurf::GetSurface(a[3]);
 458   if ( S.IsNull()) return 1;
 459
 460   Handle(Geom_Plane)   Pl = Handle(Geom_Plane)::DownCast(S);
 461   if ( Pl.IsNull()) {
 462     di << " The surface must be a plane\n";
 463     return 1;
 464   }
 465
 466   Handle(Geom_Curve) C = DrawTrSurf::GetCurve(a[2]);
 467   if ( C.IsNull()) return 1;
 468
 469   Standard_Boolean Param = Standard_True;
 470   if ((n == 5 && Draw::Atoi(a[4]) == 0) ||
 471       (n == 8 && Draw::Atoi(a[7]) == 0)) Param = Standard_False;
 472
 473   gp_Dir D;
 474
 475   if ( n == 8) {
 476     D = gp_Dir(Draw::Atof(a[4]),Draw::Atof(a[5]),Draw::Atof(a[6]));
 477   }
 478   else {
 479     D = Pl->Pln().Position().Direction();
 480   }
 481
 482   Handle(Geom_Curve) Res =
 483     GeomProjLib::ProjectOnPlane(C,Pl,D,Param);
 484
 485   DrawTrSurf::Set(a[1],Res);
 486   return 0;
 487
 488 }
 489
 490
 491 //=======================================================================
 492 //function : bisec
 493 //purpose  :
 494 //=======================================================================
 495
 496 static void solution(const Handle(GccInt_Bisec)& Bis,
 497                      const char* name,
 498                      const Standard_Integer i)
 499 {
 500   char solname[200];
 501   if ( i == 0)
 502     Sprintf(solname,"%s",name);
 503   else
 504     Sprintf(solname,"%s_%d",name,i);
 505   const char* temp = solname; // pour portage WNT
 506
 507   switch ( Bis->ArcType()) {
 508   case GccInt_Lin:
 509     DrawTrSurf::Set(temp, new Geom2d_Line(Bis->Line()));
 510     break;
 511   case GccInt_Cir:
 512     DrawTrSurf::Set(temp, new Geom2d_Circle(Bis->Circle()));
 513     break;
 514   case GccInt_Ell:
 515     DrawTrSurf::Set(temp, new Geom2d_Ellipse(Bis->Ellipse()));
 516     break;
 517   case GccInt_Par:
 518     DrawTrSurf::Set(temp, new Geom2d_Parabola(Bis->Parabola()));
 519     break;
 520   case GccInt_Hpr:
 521     DrawTrSurf::Set(temp, new Geom2d_Hyperbola(Bis->Hyperbola()));
 522     break;
 523   case GccInt_Pnt:
 524     DrawTrSurf::Set(temp, Bis->Point());
 525     break;
 526   }
 527 }
 528
 529 static Standard_Integer bisec (Draw_Interpretor& di,
 530                                Standard_Integer n, const char** a)
 531 {
 532   if (n < 4) return 1;
 533
 534   Handle(Geom2d_Curve) C1 = DrawTrSurf::GetCurve2d(a[2]);
 535   Handle(Geom2d_Curve) C2 = DrawTrSurf::GetCurve2d(a[3]);
 536   gp_Pnt2d P1,P2;
 537   Standard_Boolean ip1 = DrawTrSurf::GetPoint2d(a[2],P1);
 538   Standard_Boolean ip2 = DrawTrSurf::GetPoint2d(a[3],P2);
 539   Standard_Integer i, Compt = 0;
 540   Standard_Integer NbSol = 0;
 541
 542   if ( !C1.IsNull()) {
 543     Handle(Standard_Type) Type1 = C1->DynamicType();
 544     if ( !C2.IsNull()) {
 545       Handle(Standard_Type) Type2 = C2->DynamicType();
 546       if ( Type1 == STANDARD_TYPE(Geom2d_Line) &&
 547            Type2 == STANDARD_TYPE(Geom2d_Line)   ) {
 548         GccAna_Lin2dBisec Bis(Handle(Geom2d_Line)::DownCast(C1)->Lin2d(),
 549                               Handle(Geom2d_Line)::DownCast(C2)->Lin2d());
 550         if ( Bis.IsDone()) {
 551           char solname[200];
 552           NbSol = Bis.NbSolutions();
 553           for ( i = 1; i <= NbSol; i++) {
 554             Sprintf(solname,"%s_%d",a[1],i);
 555             const char* temp = solname; // pour portage WNT
 556             DrawTrSurf::Set(temp,new Geom2d_Line(Bis.ThisSolution(i)));
 557           }
 558         }
 559         else {
 560           di << " Bisec has failed !!\n";
 561           return 1;
 562         }
 563       }
 564       else if ( Type1 == STANDARD_TYPE(Geom2d_Line) &&
 565                 Type2 == STANDARD_TYPE(Geom2d_Circle) ) {
 566         GccAna_CircLin2dBisec
 567           Bis(Handle(Geom2d_Circle)::DownCast(C2)->Circ2d(),
 568               Handle(Geom2d_Line)::DownCast(C1)->Lin2d());
 569         if ( Bis.IsDone()) {
 570           NbSol= Bis.NbSolutions();
 571           if ( NbSol >= 2) Compt = 1;
 572           for ( i = 1; i <= NbSol; i++) {
 573             solution(Bis.ThisSolution(i),a[1],Compt);
 574             Compt++;
 575           }
 576         }
 577         else {
 578           di << " Bisec has failed !!\n";
 579           return 1;
 580         }
 581       }
 582       else if ( Type2 == STANDARD_TYPE(Geom2d_Line) &&
 583                 Type1 == STANDARD_TYPE(Geom2d_Circle) ) {
 584         GccAna_CircLin2dBisec
 585           Bis(Handle(Geom2d_Circle)::DownCast(C1)->Circ2d(),
 586               Handle(Geom2d_Line)::DownCast(C2)->Lin2d());
 587         if ( Bis.IsDone()) {
 588 //        char solname[200];
 589           NbSol = Bis.NbSolutions();
 590           if ( NbSol >= 2) Compt = 1;
 591           for ( i = 1; i <= NbSol; i++) {
 592             solution(Bis.ThisSolution(i),a[1],Compt);
 593             Compt++;
 594           }
 595         }
 596         else {
 597           di << " Bisec has failed !!\n";
 598           return 1;
 599         }
 600       }
 601       else if ( Type2 == STANDARD_TYPE(Geom2d_Circle) &&
 602                 Type1 == STANDARD_TYPE(Geom2d_Circle) ) {
 603         GccAna_Circ2dBisec
 604           Bis(Handle(Geom2d_Circle)::DownCast(C1)->Circ2d(),
 605               Handle(Geom2d_Circle)::DownCast(C2)->Circ2d());
 606         if ( Bis.IsDone()) {
 607 //        char solname[200];
 608           NbSol = Bis.NbSolutions();
 609           if ( NbSol >= 2) Compt = 1;
 610           for ( i = 1; i <= NbSol; i++) {
 611             solution(Bis.ThisSolution(i),a[1],Compt);
 612             Compt++;
 613           }
 614         }
 615         else {
 616           di << " Bisec has failed !!\n";
 617           return 1;
 618         }
 619       }
 620       else {
 621         di << " args must be line/circle/point line/circle/point\n";
 622         return 1;
 623       }
 624     }
 625     else if (ip2) {
 626       if ( Type1 == STANDARD_TYPE(Geom2d_Circle)) {
 627         GccAna_CircPnt2dBisec Bis
 628           (Handle(Geom2d_Circle)::DownCast(C1)->Circ2d(),P2);
 629         if ( Bis.IsDone()) {
 630           NbSol = Bis.NbSolutions();
 631           if ( NbSol >= 2) Compt = 1;
 632           for ( i = 1; i <= NbSol; i++) {
 633             solution(Bis.ThisSolution(i),a[1],Compt);
 634             Compt++;
 635           }
 636         }
 637         else {
 638           di << " Bisec has failed !!\n";
 639           return 1;
 640         }
 641       }
 642       else if ( Type1 == STANDARD_TYPE(Geom2d_Line)) {
 643         GccAna_LinPnt2dBisec Bis
 644           (Handle(Geom2d_Line)::DownCast(C1)->Lin2d(),P2);
 645         if ( Bis.IsDone()) {
 646           NbSol = 1;
 647           solution(Bis.ThisSolution(),a[1],0);
 648         }
 649         else {
 650           di << " Bisec has failed !!\n";
 651           return 1;
 652         }
 653       }
 654     }
 655     else {
 656       di << " the second arg must be line/circle/point \n";
 657     }
 658   }
 659   else if ( ip1) {
 660     if ( !C2.IsNull()) {
 661       Handle(Standard_Type) Type2 = C2->DynamicType();
 662       if ( Type2 == STANDARD_TYPE(Geom2d_Circle)) {
 663         GccAna_CircPnt2dBisec Bis
 664           (Handle(Geom2d_Circle)::DownCast(C2)->Circ2d(),P1);
 665         if ( Bis.IsDone()) {
 666           NbSol = Bis.NbSolutions();
 667           if ( NbSol >= 2) Compt = 1;
 668           for ( i = 1; i <= Bis.NbSolutions(); i++) {
 669             solution(Bis.ThisSolution(i),a[1],Compt);
 670             Compt++;
 671           }
 672         }
 673         else {
 674           di << " Bisec has failed !!\n";
 675           return 1;
 676         }
 677       }
 678       else if ( Type2 == STANDARD_TYPE(Geom2d_Line)) {
 679         GccAna_LinPnt2dBisec Bis
 680           (Handle(Geom2d_Line)::DownCast(C2)->Lin2d(),P1);
 681         if ( Bis.IsDone()) {
 682           NbSol = 1;
 683           solution(Bis.ThisSolution(),a[1],0);
 684         }
 685         else {
 686           di << " Bisec has failed !!\n";
 687           return 1;
 688         }
 689       }
 690     }
 691     else if (ip2) {
 692       GccAna_Pnt2dBisec Bis(P1,P2);
 693       if ( Bis.HasSolution()) {
 694         NbSol = 1;
 695         DrawTrSurf::Set(a[1],new Geom2d_Line(Bis.ThisSolution()));
 696       }
 697       else {
 698         di << " Bisec has failed !!\n";
 699         return 1;
 700       }
 701     }
 702     else {
 703       di << " the second arg must be line/circle/point \n";
 704       return 1;
 705     }
 706   }
 707   else {
 708     di << " args must be line/circle/point line/circle/point\n";
 709     return 1;
 710   }
 711
 712   if ( NbSol >= 2) {
 713     di << "There are " << NbSol << " Solutions.\n";
 714   }
 715   else {
 716     di << "There is " << NbSol << " Solution.\n";
 717   }
 718
 719   return 0;
 720 }
 721
 722 //=======================================================================
 723 //function : cmovetangent
 724 //purpose  :
 725 //=======================================================================
 726
 727 static Standard_Integer movelaw (Draw_Interpretor& di, Standard_Integer n, const char** a)
 728 {
 729   Standard_Integer
 730     ii,
 731     condition=0,
 732     error_status ;
 733   Standard_Real u,
 734     x,
 735     tolerance,
 736     tx ;
 737
 738   u = Draw::Atof(a[2]);
 739   x = Draw::Atof(a[3]);
 740   tolerance = 1.0e-5 ;
 741   if (n < 5) {
 742     return 1 ;
 743   }
 744   Handle(Geom2d_BSplineCurve) G2 = DrawTrSurf::GetBSplineCurve2d(a[1]);
 745   if (!G2.IsNull()) {
 746     tx = Draw::Atof(a[4]) ;
 747     if (n == 6) {
 748       condition = Max(Draw::Atoi(a[5]), -1)  ;
 749       condition = Min(condition, G2->Degree()-1) ;
 750     }
 751     TColgp_Array1OfPnt2d   curve_poles(1,G2->NbPoles()) ;
 752     TColStd_Array1OfReal    law_poles(1,G2->NbPoles()) ;
 753     TColStd_Array1OfReal    law_knots(1,G2->NbKnots()) ;
 754     TColStd_Array1OfInteger law_mults(1,G2->NbKnots()) ;
 755
 756     G2->Knots(law_knots) ;
 757     G2->Multiplicities(law_mults) ;
 758     G2->Poles(curve_poles) ;
 759     for (ii = 1 ; ii <= G2->NbPoles() ; ii++) {
 760       law_poles(ii) = curve_poles(ii).Coord(2) ;
 761     }
 762
 763     Law_BSpline  a_law(law_poles,
 764       law_knots,
 765       law_mults,
 766       G2->Degree(),
 767       Standard_False) ;
 768
 769     a_law.MovePointAndTangent(u,
 770       x,
 771       tx,
 772       tolerance,
 773       condition,
 774       condition,
 775       error_status) ;
 776
 777     for (ii = 1 ; ii <= G2->NbPoles() ; ii++) {
 778       curve_poles(ii).SetCoord(2,a_law.Pole(ii)) ;
 779       G2->SetPole(ii,curve_poles(ii)) ;
 780     }
 781
 782
 783     if (! error_status) {
 784       Draw::Repaint();
 785     }
 786     else {
 787       di << "Not enough degree of freedom increase degree please\n";
 788     }
 789
 790
 791   }
 792   return 0;
 793 }
 794
 795
 796 //Static method computing deviation of curve and polyline
 797 #include <math_PSO.hxx>
 798 #include <math_PSOParticlesPool.hxx>
 799 #include <math_MultipleVarFunction.hxx>
 800 #include <math_BrentMinimum.hxx>
 801
 802 static Standard_Real CompLocalDev(const Adaptor3d_Curve& theCurve,
 803                                   const Standard_Real u1, const Standard_Real u2);
 804
 805 static void ComputeDeviation(const Adaptor3d_Curve& theCurve,
 806                              const Handle(Geom_BSplineCurve)& thePnts,
 807                              Standard_Real& theDmax,
 808                              Standard_Real& theUfMax,
 809                              Standard_Real& theUlMax,
 810                              Standard_Integer& theImax)
 811 {
 812   theDmax = 0.;
 813   theUfMax = 0.;
 814   theUlMax = 0.;
 815   theImax = 0;
 816
 817   //take knots
 818   Standard_Integer nbp = thePnts->NbKnots();
 819   TColStd_Array1OfReal aKnots(1, nbp);
 820   thePnts->Knots(aKnots);
 821
 822   Standard_Integer i;
 823   for(i = 1; i < nbp; ++i)
 824   {
 825     Standard_Real u1 = aKnots(i), u2 = aKnots(i+1);
 826     Standard_Real d = CompLocalDev(theCurve, u1, u2);
 827     if(d > theDmax)
 828     {
 829       theDmax = d;
 830       theImax = i;
 831       theUfMax = u1;
 832       theUlMax = u2;
 833     }
 834   }
 835 }
 836
 837 Standard_Real CompLocalDev(const Adaptor3d_Curve& theCurve,
 838                            const Standard_Real u1, const Standard_Real u2)
 839 {
 840   math_Vector aLowBorder(1,1);
 841   math_Vector aUppBorder(1,1);
 842   math_Vector aSteps(1,1);
 843   //
 844   aLowBorder(1) = u1;
 845   aUppBorder(1) = u2;
 846   aSteps(1) =(aUppBorder(1) - aLowBorder(1)) * 0.01; // Run PSO on even distribution with 100 points.
 847   //
 848   GCPnts_DistFunction aFunc1(theCurve,  u1, u2);
 849   //
 850   Standard_Real aValue;
 851   math_Vector aT(1,1);
 852   GCPnts_DistFunctionMV aFunc(aFunc1);
 853
 854   math_PSO aFinder(&aFunc, aLowBorder, aUppBorder, aSteps); // Choose 32 best points from 100 above.
 855   aFinder.Perform(aSteps, aValue, aT);
 856   Standard_Real d = 0.;
 857
 858   Standard_Real d1, d2;
 859   Standard_Real x1 = Max(u1, aT(1) - aSteps(1));
 860   Standard_Boolean Ok = aFunc1.Value(x1, d1);
 861   if(!Ok)
 862   {
 863     return Sqrt(-aValue);
 864   }
 865   Standard_Real x2 = Min(u2, aT(1) + aSteps(1));
 866   Ok = aFunc1.Value(x2, d2);
 867   if(!Ok)
 868   {
 869     return Sqrt(-aValue);
 870   }
 871   if(!(d1 > aValue && d2 > aValue))
 872   {
 873     Standard_Real dmin = Min(d1, Min(aValue, d2));
 874     return Sqrt(-dmin);
 875   }
 876
 877   math_BrentMinimum anOptLoc(Precision::PConfusion());
 878   anOptLoc.Perform(aFunc1, x1, aT(1), x2);
 879
 880   if (anOptLoc.IsDone())
 881   {
 882     d = -anOptLoc.Minimum();
 883   }
 884   else
 885   {
 886     d = -aValue;
 887   }
 888   return Sqrt(d);
 889 }
 890
 891 //=======================================================================
 892 //function : crvpoints
 893 //purpose  :
 894 //=======================================================================
 895
 896 static Standard_Integer crvpoints (Draw_Interpretor& di, Standard_Integer /*n*/, const char** a)
 897 {
 898   Standard_Integer i, nbp;
 899   Standard_Real defl;
 900
 901   Handle(Adaptor3d_Curve) aHCurve;
 902   Handle(Geom_Curve) C = DrawTrSurf::GetCurve(a[2]);
 903   if (C.IsNull())
 904   {
 905     // try getting a wire
 906     TopoDS_Wire aWire = TopoDS::Wire(DBRep::Get(a[2], TopAbs_WIRE));
 907     if (aWire.IsNull())
 908     {
 909       Message::SendFail() << "cannot evaluate the argument " << a[2] << " as a curve";
 910       return 1;
 911     }
 912     BRepAdaptor_CompCurve aCompCurve(aWire);
 913     aHCurve = new BRepAdaptor_CompCurve(aCompCurve);
 914   }
 915   else
 916   {
 917     aHCurve = new GeomAdaptor_Curve(C);
 918   }
 919
 920   defl = Draw::Atof(a[3]);
 921
 922   GCPnts_QuasiUniformDeflection PntGen (*aHCurve, defl);
 923
 924   if(!PntGen.IsDone()) {
 925     di << "Points generation failed\n";
 926     return 1;
 927   }
 928
 929   nbp = PntGen.NbPoints();
 930   di << "Nb points : " << nbp << "\n";
 931
 932   TColgp_Array1OfPnt aPoles(1, nbp);
 933   TColStd_Array1OfReal aKnots(1, nbp);
 934   TColStd_Array1OfInteger aMults(1, nbp);
 935
 936   for(i = 1; i <= nbp; ++i) {
 937     aPoles(i) = PntGen.Value(i);
 938     aKnots(i) = PntGen.Parameter(i);
 939     aMults(i) = 1;
 940   }
 941
 942   aMults(1) = 2;
 943   aMults(nbp) = 2;
 944
 945   Handle(Geom_BSplineCurve) aPnts = new Geom_BSplineCurve(aPoles, aKnots, aMults, 1);
 946   Handle(DrawTrSurf_BSplineCurve) aDrCrv = new DrawTrSurf_BSplineCurve(aPnts);
 947
 948   aDrCrv->ClearPoles();
 949   Draw_Color aKnColor(Draw_or);
 950   aDrCrv->SetKnotsColor(aKnColor);
 951   aDrCrv->SetKnotsShape(Draw_Plus);
 952
 953   Draw::Set(a[1], aDrCrv);
 954
 955   Standard_Real dmax = 0., ufmax = 0., ulmax = 0.;
 956   Standard_Integer imax = 0;
 957
 958   //check deviation
 959   ComputeDeviation (*aHCurve, aPnts, dmax, ufmax, ulmax, imax);
 960   di << "Max defl: " << dmax << " " << ufmax << " " << ulmax << " " << imax << "\n";
 961
 962   return 0;
 963 }
 964
 965 //=======================================================================
 966 //function : crvtpoints
 967 //purpose  :
 968 //=======================================================================
 969
 970 static Standard_Integer crvtpoints (Draw_Interpretor& di, Standard_Integer n, const char** a)
 971 {
 972   Standard_Integer i, nbp, aMinPntsNb = 2;
 973   Standard_Real defl, angle = Precision::Angular();
 974
 975   Handle(Adaptor3d_Curve) aHCurve;
 976   Handle(Geom_Curve) C = DrawTrSurf::GetCurve(a[2]);
 977   if (C.IsNull())
 978   {
 979     // try getting a wire
 980     TopoDS_Wire aWire = TopoDS::Wire(DBRep::Get(a[2], TopAbs_WIRE));
 981     if (aWire.IsNull())
 982     {
 983       Message::SendFail() << "cannot evaluate the argument " << a[2] << " as a curve";
 984       return 1;
 985     }
 986     BRepAdaptor_CompCurve aCompCurve(aWire);
 987     aHCurve = new BRepAdaptor_CompCurve(aCompCurve);
 988   }
 989   else
 990   {
 991     aHCurve = new GeomAdaptor_Curve(C);
 992   }
 993   defl = Draw::Atof(a[3]);
 994
 995   if(n > 4)
 996     angle = Draw::Atof(a[4]);
 997
 998   if(n > 5)
 999     aMinPntsNb = Draw::Atoi (a[5]);
1000
1001   GCPnts_TangentialDeflection PntGen (*aHCurve, angle, defl, aMinPntsNb);
1002
1003   nbp = PntGen.NbPoints();
1004   di << "Nb points : " << nbp << "\n";
1005
1006   TColgp_Array1OfPnt aPoles(1, nbp);
1007   TColStd_Array1OfReal aKnots(1, nbp);
1008   TColStd_Array1OfInteger aMults(1, nbp);
1009
1010   for(i = 1; i <= nbp; ++i) {
1011     aPoles(i) = PntGen.Value(i);
1012     aKnots(i) = PntGen.Parameter(i);
1013     aMults(i) = 1;
1014   }
1015
1016   aMults(1) = 2;
1017   aMults(nbp) = 2;
1018
1019   Handle(Geom_BSplineCurve) aPnts = new Geom_BSplineCurve(aPoles, aKnots, aMults, 1);
1020   Handle(DrawTrSurf_BSplineCurve) aDrCrv = new DrawTrSurf_BSplineCurve(aPnts);
1021
1022   aDrCrv->ClearPoles();
1023   Draw_Color aKnColor(Draw_or);
1024   aDrCrv->SetKnotsColor(aKnColor);
1025   aDrCrv->SetKnotsShape(Draw_Plus);
1026
1027   Draw::Set(a[1], aDrCrv);
1028
1029   Standard_Real dmax = 0., ufmax = 0., ulmax = 0.;
1030   Standard_Integer imax = 0;
1031
1032   //check deviation
1033   ComputeDeviation (*aHCurve, aPnts, dmax, ufmax, ulmax, imax);
1034   //
1035   di << "Max defl: " << dmax << " " << ufmax << " " << ulmax << " " << imax << "\n";
1036
1037   return 0;
1038 }
1039 //=======================================================================
1040 //function : uniformAbscissa
1041 //purpose  : epa test (TATA-06-002 (Problem with GCPnts_UniformAbscissa class)
1042 //=======================================================================
1043 static Standard_Integer uniformAbscissa (Draw_Interpretor& di, Standard_Integer n, const char** a)
1044 {
1045   if( n != 3 )
1046     return 1;
1047
1048   /*Handle(Geom_BSplineCurve) ellip;
1049   ellip = DrawTrSurf::GetBSplineCurve(a[1]);
1050   if (ellip.IsNull())
1051   {
1052     di << " BSpline is NULL  \n";
1053     return 1;
1054   }*/
1055
1056   Handle(Geom_Curve) ellip;
1057   ellip = DrawTrSurf::GetCurve(a[1]);
1058   if (ellip.IsNull())
1059   {
1060     di << " Curve is NULL  \n";
1061     return 1;
1062   }
1063
1064   Standard_Integer nocp;
1065   nocp = Draw::Atoi(a[2]);
1066   if(nocp < 2)
1067     return 1;
1068
1069
1070   //test nbPoints for Geom_Ellipse
1071
1072   try
1073   {
1074     GeomLProp_CLProps Prop(ellip,2,Precision::Intersection());
1075     Prop.SetCurve(ellip);
1076
1077     GeomAdaptor_Curve GAC(ellip);
1078     di<<"Type Of curve: "<<GAC.GetType()<<"\n";
1079     Standard_Real Tol = Precision::Confusion();
1080     Standard_Real L;
1081
1082     L = GCPnts_AbscissaPoint::Length(GAC, GAC.FirstParameter(), GAC.LastParameter(), Tol);
1083     di<<"Ellipse length = "<<L<<"\n";
1084     Standard_Real Abscissa = L/(nocp-1);
1085     di << " CUR : Abscissa " << Abscissa << "\n";
1086
1087     GCPnts_UniformAbscissa myAlgo(GAC, Abscissa, ellip->FirstParameter(), ellip->LastParameter());
1088     if ( myAlgo.IsDone() )
1089     {
1090       di << " CasCurve  - nbpoints " << myAlgo.NbPoints() << "\n";
1091       for(Standard_Integer i = 1; i<= myAlgo.NbPoints(); i++ )
1092         di << i <<" points = " << myAlgo.Parameter( i ) << "\n";
1093     }
1094   }
1095
1096   catch (Standard_Failure const&)
1097   {
1098     di << " Standard Failure  \n";
1099   }
1100   return 0;
1101 }
1102
1103 //=======================================================================
1104 //function : EllipsUniformAbscissa
1105 //purpose  : epa test (TATA-06-002 (Problem with GCPnts_UniformAbscissa class)
1106 //=======================================================================
1107 static Standard_Integer EllipsUniformAbscissa (Draw_Interpretor& di, Standard_Integer n, const char** a)
1108 {
1109   if( n != 4 )
1110     return 1;
1111
1112   Standard_Real R1;
1113   R1 = Draw::Atof(a[1]);
1114   Standard_Real R2;
1115   R2 = Draw::Atof(a[2]);
1116
1117   Standard_Integer nocp;
1118   nocp = Draw::Atoi(a[3]);
1119   if(nocp < 2)
1120     return 1;
1121
1122   //test nbPoints for Geom_Ellipse
1123   Handle(Geom_Ellipse) ellip;
1124
1125
1126   try
1127   {
1128     gp_Pnt location;
1129     location = gp_Pnt( 0.0, 0.0, 0.0);
1130     gp_Dir main_direction(0.0, 0.0, 1.0);
1131
1132     gp_Dir x_direction(1.0, 0.0, 0.0);
1133     gp_Ax2 mainaxis( location, main_direction);
1134
1135     mainaxis.SetXDirection(x_direction);
1136     ellip = new Geom_Ellipse(mainaxis,R1, R2);
1137
1138     BRepBuilderAPI_MakeEdge curve_edge(ellip);
1139     TopoDS_Edge edge_curve = curve_edge.Edge();
1140
1141     DBRep::Set("Ellipse",edge_curve);
1142   }
1143
1144   catch(Standard_Failure const&)
1145   {
1146     di << " Standard Failure  \n";
1147   }
1148
1149   try
1150   {
1151     GeomLProp_CLProps Prop(ellip,2,Precision::Intersection());
1152     Prop.SetCurve(ellip);
1153
1154     GeomAdaptor_Curve GAC(ellip);
1155     di<<"Type Of curve: "<<GAC.GetType()<<"\n";
1156     Standard_Real Tol = Precision::Confusion();
1157     Standard_Real L;
1158
1159     L = GCPnts_AbscissaPoint::Length(GAC, GAC.FirstParameter(), GAC.LastParameter(), Tol);
1160     di<<"Ellipse length = "<<L<<"\n";
1161     Standard_Real Abscissa = L/(nocp-1);
1162     di << " CUR : Abscissa " << Abscissa << "\n";
1163
1164     GCPnts_UniformAbscissa myAlgo(GAC, Abscissa, ellip->FirstParameter(), ellip->LastParameter());
1165     if ( myAlgo.IsDone() )
1166     {
1167       di << " CasCurve  - nbpoints " << myAlgo.NbPoints() << "\n";
1168       for(Standard_Integer i = 1; i<= myAlgo.NbPoints(); i++ )
1169         di << i <<" points = " << myAlgo.Parameter( i ) << "\n";
1170     }
1171   }
1172
1173   catch (Standard_Failure const&)
1174   {
1175     di << " Standard Failure  \n";
1176   }
1177   return 0;
1178 }
1179
1180 //=======================================================================
1181 //function : discrCurve
1182 //purpose  :
1183 //=======================================================================
1184 static Standard_Integer discrCurve(Draw_Interpretor& di, Standard_Integer theArgNb, const char** theArgVec)
1185 {
1186   if (theArgNb < 3)
1187   {
1188     di << "Invalid number of parameters.\n";
1189     return 1;
1190   }
1191
1192   Handle(Geom_Curve) aCurve = DrawTrSurf::GetCurve(theArgVec[2]);
1193   if (aCurve.IsNull())
1194   {
1195     di << "Curve is NULL.\n";
1196     return 1;
1197   }
1198
1199   Standard_Integer aSrcNbPnts = 0;
1200   Standard_Boolean isUniform = Standard_False;
1201   for (Standard_Integer anArgIter = 3; anArgIter < theArgNb; ++anArgIter)
1202   {
1203     TCollection_AsciiString anArg     (theArgVec[anArgIter]);
1204     TCollection_AsciiString anArgCase (anArg);
1205     anArgCase.LowerCase();
1206     if (anArgCase == "nbpnts")
1207     {
1208       if (++anArgIter >= theArgNb)
1209       {
1210         di << "Value for argument '" << anArg << "' is absent.\n";
1211         return 1;
1212       }
1213
1214       aSrcNbPnts = Draw::Atoi (theArgVec[anArgIter]);
1215     }
1216     else if (anArgCase == "uniform")
1217     {
1218       if (++anArgIter >= theArgNb)
1219       {
1220         di << "Value for argument '" << anArg << "' is absent.\n";
1221         return 1;
1222       }
1223
1224       isUniform = (Draw::Atoi (theArgVec[anArgIter]) == 1);
1225     }
1226     else
1227     {
1228       di << "Invalid argument '" << anArg << "'.\n";
1229       return 1;
1230     }
1231   }
1232
1233   if (aSrcNbPnts < 2)
1234   {
1235     di << "Invalid count of points.\n";
1236     return 1;
1237   }
1238
1239   if (!isUniform)
1240   {
1241     di << "Invalid type of discretization.\n";
1242     return 1;
1243   }
1244
1245   GeomAdaptor_Curve aCurveAdaptor(aCurve);
1246   GCPnts_UniformAbscissa aSplitter(aCurveAdaptor, aSrcNbPnts, Precision::Confusion());
1247   if (!aSplitter.IsDone())
1248   {
1249     di << "Error: Invalid result.\n";
1250     return 0;
1251   }
1252
1253   const Standard_Integer aDstNbPnts = aSplitter.NbPoints();
1254
1255   if (aDstNbPnts < 2)
1256   {
1257     di << "Error: Invalid result.\n";
1258     return 0;
1259   }
1260
1261   TColgp_Array1OfPnt aPoles(1, aDstNbPnts);
1262   TColStd_Array1OfReal aKnots(1, aDstNbPnts);
1263   TColStd_Array1OfInteger aMultiplicities(1, aDstNbPnts);
1264
1265   for (Standard_Integer aPntIter = 1; aPntIter <= aDstNbPnts; ++aPntIter)
1266   {
1267     aPoles.ChangeValue(aPntIter) = aCurveAdaptor.Value(aSplitter.Parameter(aPntIter));
1268     aKnots.ChangeValue(aPntIter) = (aPntIter - 1) / (aDstNbPnts - 1.0);
1269     aMultiplicities.ChangeValue(aPntIter) = 1;
1270   }
1271   aMultiplicities.ChangeValue(1) = 2;
1272   aMultiplicities.ChangeValue(aDstNbPnts) = 2;
1273
1274   Handle(Geom_BSplineCurve) aPolyline =
1275     new Geom_BSplineCurve(aPoles, aKnots, aMultiplicities, 1);
1276   DrawTrSurf::Set(theArgVec[1], aPolyline);
1277
1278   return 0;
1279 }
1280
1281 //=======================================================================
1282 //function : mypoints
1283 //purpose  :
1284 //=======================================================================
1285
1286 static Standard_Integer mypoints (Draw_Interpretor& di, Standard_Integer /*n*/, const char** a)
1287 {
1288   Standard_Integer i, nbp;
1289   Standard_Real defl;
1290
1291   Handle(Geom_Curve) C = DrawTrSurf::GetCurve(a[2]);
1292   defl = Draw::Atof(a[3]);
1293   Handle(Geom_BSplineCurve) aBS (Handle(Geom_BSplineCurve)::DownCast(C));
1294
1295   if(aBS.IsNull()) return 1;
1296
1297   Standard_Integer ui1 = aBS->FirstUKnotIndex();
1298   Standard_Integer ui2 = aBS->LastUKnotIndex();
1299
1300   Standard_Integer nbsu = ui2-ui1+1; nbsu += (nbsu - 1) * (aBS->Degree()-1);
1301
1302   TColStd_Array1OfReal anUPars(1, nbsu);
1303   TColStd_Array1OfBoolean anUFlg(1, nbsu);
1304
1305   Standard_Integer j, k, nbi;
1306   Standard_Real t1, t2, dt;
1307
1308   //Filling of sample parameters
1309   nbi = aBS->Degree();
1310   k = 0;
1311   t1 = aBS->Knot(ui1);
1312   for(i = ui1+1; i <= ui2; ++i) {
1313     t2 = aBS->Knot(i);
1314     dt = (t2 - t1)/nbi;
1315     j = 1;
1316     do {
1317       ++k;
1318       anUPars(k) = t1;
1319       anUFlg(k) = Standard_False;
1320       t1 += dt;
1321     }
1322     while (++j <= nbi);
1323     t1 = t2;
1324   }
1325   ++k;
1326   anUPars(k) = t1;
1327
1328   Standard_Integer l;
1329   defl *= defl;
1330
1331   j = 1;
1332   anUFlg(1) = Standard_True;
1333   anUFlg(nbsu) = Standard_True;
1334   Standard_Boolean bCont = Standard_True;
1335   while (j < nbsu-1 && bCont) {
1336     t2 = anUPars(j);
1337     gp_Pnt p1 = aBS->Value(t2);
1338     for(k = j+2; k <= nbsu; ++k) {
1339       t2 = anUPars(k);
1340       gp_Pnt p2 = aBS->Value(t2);
1341       gce_MakeLin MkLin(p1, p2);
1342       const gp_Lin& lin = MkLin.Value();
1343       Standard_Boolean ok = Standard_True;
1344       for(l = j+1; l < k; ++l) {
1345         if(anUFlg(l)) continue;
1346         gp_Pnt pp =  aBS->Value(anUPars(l));
1347         Standard_Real d = lin.SquareDistance(pp);
1348
1349         if(d <= defl) continue;
1350
1351         ok = Standard_False;
1352         break;
1353       }
1354
1355
1356       if(!ok) {
1357         j = k - 1;
1358         anUFlg(j) = Standard_True;
1359         break;
1360       }
1361
1362     }
1363
1364     if(k >= nbsu) bCont = Standard_False;
1365   }
1366
1367   nbp = 0;
1368   for(i = 1; i <= nbsu; ++i) {
1369     if(anUFlg(i)) nbp++;
1370   }
1371
1372   TColgp_Array1OfPnt aPoles(1, nbp);
1373   TColStd_Array1OfReal aKnots(1, nbp);
1374   TColStd_Array1OfInteger aMults(1, nbp);
1375   j = 0;
1376   for(i = 1; i <= nbsu; ++i) {
1377     if(anUFlg(i)) {
1378       ++j;
1379       aKnots(j) = anUPars(i);
1380       aMults(j) = 1;
1381       aPoles(j) = aBS->Value(aKnots(j));
1382     }
1383   }
1384
1385   aMults(1) = 2;
1386   aMults(nbp) = 2;
1387
1388   Handle(Geom_BSplineCurve) aPnts = new Geom_BSplineCurve(aPoles, aKnots, aMults, 1);
1389   Handle(DrawTrSurf_BSplineCurve) aDrCrv = new DrawTrSurf_BSplineCurve(aPnts);
1390
1391   aDrCrv->ClearPoles();
1392   Draw_Color aKnColor(Draw_or);
1393   aDrCrv->SetKnotsColor(aKnColor);
1394   aDrCrv->SetKnotsShape(Draw_Plus);
1395
1396   Draw::Set(a[1], aDrCrv);
1397
1398   Standard_Real dmax = 0., ufmax = 0., ulmax = 0.;
1399   Standard_Integer imax = 0;
1400
1401   ComputeDeviation(GeomAdaptor_Curve(C),aPnts,dmax,ufmax,ulmax,imax);
1402   di << "Max defl: " << dmax << " " << ufmax << " " << ulmax << " " << imax << "\n";
1403
1404   return 0;
1405 }
1406
1407
1408
1409 //=======================================================================
1410 //function : surfpoints
1411 //purpose  :
1412 //=======================================================================
1413
1414 static Standard_Integer surfpoints (Draw_Interpretor& /*di*/, Standard_Integer /*n*/, const char** a)
1415 {
1416   Standard_Integer i;
1417   Standard_Real defl;
1418
1419   Handle(Geom_Surface) S = DrawTrSurf::GetSurface(a[2]);
1420   defl = Draw::Atof(a[3]);
1421
1422   Handle(GeomAdaptor_Surface) AS = new GeomAdaptor_Surface(S);
1423
1424   Handle(Adaptor3d_TopolTool) aTopTool = new Adaptor3d_TopolTool(AS);
1425
1426   aTopTool->SamplePnts(defl, 10, 10);
1427
1428   Standard_Integer nbpu = aTopTool->NbSamplesU();
1429   Standard_Integer nbpv = aTopTool->NbSamplesV();
1430   TColStd_Array1OfReal Upars(1, nbpu), Vpars(1, nbpv);
1431   aTopTool->UParameters(Upars);
1432   aTopTool->VParameters(Vpars);
1433
1434   TColgp_Array2OfPnt aPoles(1, nbpu, 1, nbpv);
1435   TColStd_Array1OfReal anUKnots(1, nbpu);
1436   TColStd_Array1OfReal aVKnots(1, nbpv);
1437   TColStd_Array1OfInteger anUMults(1, nbpu);
1438   TColStd_Array1OfInteger aVMults(1, nbpv);
1439
1440   Standard_Integer j;
1441   for(i = 1; i <= nbpu; ++i) {
1442     anUKnots(i) = Upars(i);
1443     anUMults(i) = 1;
1444     for(j = 1; j <= nbpv; ++j) {
1445       aVKnots(j) = Vpars(j);
1446       aVMults(j) = 1;
1447       aPoles(i,j) = S->Value(anUKnots(i),aVKnots(j));
1448     }
1449   }
1450
1451   anUMults(1) = 2;
1452   anUMults(nbpu) = 2;
1453   aVMults(1) = 2;
1454   aVMults(nbpv) = 2;
1455
1456   Handle(Geom_BSplineSurface) aPnts = new Geom_BSplineSurface(aPoles, anUKnots,  aVKnots,
1457                                                               anUMults, aVMults, 1, 1);
1458   Handle(DrawTrSurf_BSplineSurface) aDrSurf = new DrawTrSurf_BSplineSurface(aPnts);
1459
1460   aDrSurf->ClearPoles();
1461   Draw_Color aKnColor(Draw_or);
1462   aDrSurf->SetKnotsColor(aKnColor);
1463   aDrSurf->SetKnotsShape(Draw_Plus);
1464
1465   Draw::Set(a[1], aDrSurf);
1466
1467
1468   return 0;
1469 }
1470
1471
1472
1473 //=======================================================================
1474 //function : intersect
1475 //purpose  :
1476 //=======================================================================
1477 static Standard_Integer intersection (Draw_Interpretor& di,
1478                                       Standard_Integer n, const char** a)
1479 {
1480   if (n < 4)
1481     return 1;
1482
1483   //
1484   Handle(Geom_Curve) GC1;
1485   Handle(Geom_Surface) GS1 = DrawTrSurf::GetSurface(a[2]);
1486   if (GS1.IsNull())
1487   {
1488     GC1 = DrawTrSurf::GetCurve(a[2]);
1489     if (GC1.IsNull())
1490       return 1;
1491   }
1492
1493   //
1494   Handle(Geom_Surface) GS2 = DrawTrSurf::GetSurface(a[3]);
1495   if (GS2.IsNull())
1496     return 1;
1497
1498   //
1499   Standard_Real tol = Precision::Confusion();
1500   if (n == 5 || n == 9 || n == 13 || n == 17)
1501     tol = Draw::Atof(a[n-1]);
1502
1503   //
1504   Handle(Geom_Curve) Result;
1505   gp_Pnt             Point;
1506
1507   //
1508   if (GC1.IsNull())
1509   {
1510     GeomInt_IntSS Inters;
1511     //
1512     // Surface Surface
1513     if (n <= 5)
1514     {
1515       // General case
1516       Inters.Perform(GS1,GS2,tol,Standard_True);
1517     }
1518     else if (n == 8 || n == 9 || n == 12 || n == 13 || n == 16 || n == 17)
1519     {
1520       Standard_Boolean useStart = Standard_True, useBnd = Standard_True;
1521       Standard_Integer ista1=0,ista2=0,ibnd1=0,ibnd2=0;
1522       Standard_Real UVsta[4];
1523       Handle(GeomAdaptor_Surface) AS1,AS2;
1524
1525       //
1526       if (n <= 9)          // user starting point
1527       {
1528         useBnd = Standard_False;
1529         ista1 = 4;
1530         ista2 = 7;
1531       }
1532       else if (n <= 13)   // user bounding
1533       {
1534         useStart = Standard_False;
1535         ibnd1 = 4; ibnd2 = 11;
1536       }
1537       else        // both user starting point and bounding
1538       {
1539         ista1 = 4; ista2 = 7;
1540         ibnd1 = 8; ibnd2 = 15;
1541       }
1542
1543       if (useStart)
1544       {
1545         for (Standard_Integer i=ista1; i <= ista2; i++)
1546         {
1547           UVsta[i-ista1] = Draw::Atof(a[i]);
1548         }
1549       }
1550
1551       if (useBnd)
1552       {
1553         Standard_Real UVbnd[8];
1554         for (Standard_Integer i=ibnd1; i <= ibnd2; i++)
1555           UVbnd[i-ibnd1] = Draw::Atof(a[i]);
1556
1557         AS1 = new GeomAdaptor_Surface(GS1,UVbnd[0],UVbnd[1],UVbnd[2],UVbnd[3]);
1558         AS2 = new GeomAdaptor_Surface(GS2,UVbnd[4],UVbnd[5],UVbnd[6],UVbnd[7]);
1559       }
1560
1561       //
1562       if (useStart && !useBnd)
1563       {
1564         Inters.Perform(GS1,GS2,tol,UVsta[0],UVsta[1],UVsta[2],UVsta[3]);
1565       }
1566       else if (!useStart && useBnd)
1567       {
1568         Inters.Perform(AS1,AS2,tol);
1569       }
1570       else
1571       {
1572         Inters.Perform(AS1,AS2,tol,UVsta[0],UVsta[1],UVsta[2],UVsta[3]);
1573       }
1574     }//else if (n == 8 || n == 9 || n == 12 || n == 13 || n == 16 || n == 17)
1575     else
1576     {
1577       di<<"incorrect number of arguments\n";
1578       return 1;
1579     }
1580
1581     //
1582     if (!Inters.IsDone())
1583     {
1584       di<<"No intersections found!\n";
1585
1586       return 1;
1587     }
1588
1589     //
1590     char buf[1024];
1591     Standard_Integer i, aNbLines, aNbPoints;
1592
1593     //
1594     aNbLines = Inters.NbLines();
1595     if (aNbLines >= 2)
1596     {
1597       for (i=1; i<=aNbLines; ++i)
1598       {
1599         Sprintf(buf, "%s_%d",a[1],i);
1600         di << buf << " ";
1601         Result = Inters.Line(i);
1602         const char* temp = buf;
1603         DrawTrSurf::Set(temp,Result);
1604       }
1605     }
1606     else if (aNbLines == 1)
1607     {
1608       Result = Inters.Line(1);
1609       Sprintf(buf,"%s",a[1]);
1610       di << buf << " ";
1611       DrawTrSurf::Set(a[1],Result);
1612     }
1613
1614     //
1615     aNbPoints=Inters.NbPoints();
1616     for (i=1; i<=aNbPoints; ++i)
1617     {
1618       Point=Inters.Point(i);
1619       Sprintf(buf,"%s_p_%d",a[1],i);
1620       di << buf << " ";
1621       const char* temp = buf;
1622       DrawTrSurf::Set(temp, Point);
1623     }
1624   }// if (GC1.IsNull())
1625   else
1626   {
1627     // Curve Surface
1628     GeomAPI_IntCS Inters(GC1,GS2);
1629
1630     //
1631     if (!Inters.IsDone())
1632     {
1633       di<<"No intersections found!\n";
1634       return 1;
1635     }
1636
1637     Standard_Integer nblines = Inters.NbSegments();
1638     Standard_Integer nbpoints = Inters.NbPoints();
1639
1640     char newname[1024];
1641
1642     if ( (nblines+nbpoints) >= 2)
1643     {
1644       Standard_Integer i;
1645       Standard_Integer Compt = 1;
1646
1647       if(nblines >= 1)
1648         std::cout << "   Lines: " << std::endl;
1649
1650       for (i = 1; i <= nblines; i++, Compt++)
1651       {
1652         Sprintf(newname,"%s_%d",a[1],Compt);
1653         di << newname << " ";
1654         Result = Inters.Segment(i);
1655         const char* temp = newname; // pour portage WNT
1656         DrawTrSurf::Set(temp,Result);
1657       }
1658
1659       if(nbpoints >= 1)
1660         std::cout << "   Points: " << std::endl;
1661
1662       const Standard_Integer imax = nblines+nbpoints;
1663
1664       for (/*i = 1*/; i <= imax; i++, Compt++)
1665       {
1666         Sprintf(newname,"%s_%d",a[1],i);
1667         di << newname << " ";
1668         Point = Inters.Point(i);
1669         const char* temp = newname; // pour portage WNT
1670         DrawTrSurf::Set(temp,Point);
1671       }
1672     }
1673     else if (nblines == 1)
1674     {
1675       Result = Inters.Segment(1);
1676       Sprintf(newname,"%s",a[1]);
1677       di << newname << " ";
1678       DrawTrSurf::Set(a[1],Result);
1679     }
1680     else if (nbpoints == 1)
1681     {
1682       Point = Inters.Point(1);
1683       Sprintf(newname,"%s",a[1]);
1684       di << newname << " ";
1685       DrawTrSurf::Set(a[1],Point);
1686     }
1687   }
1688
1689   dout.Flush();
1690   return 0;
1691 }
1692
1693 //=======================================================================
1694 //function : GetCurveContinuity
1695 //purpose  : Returns the continuity of the given curve
1696 //=======================================================================
1697 static Standard_Integer GetCurveContinuity( Draw_Interpretor& theDI,
1698                                             Standard_Integer theNArg,
1699                                             const char** theArgv)
1700 {
1701   if(theNArg != 2)
1702   {
1703     theDI << "Use: getcurvcontinuity {curve or 2dcurve} \n";
1704     return 1;
1705   }
1706
1707   char aContName[7][3] = {"C0",   //0
1708                           "G1",   //1
1709                           "C1",   //2
1710                           "G2",   //3
1711                           "C2",   //4
1712                           "C3",   //5
1713                           "CN"};  //6
1714
1715   Handle(Geom2d_Curve) GC2d;
1716   Handle(Geom_Curve) GC3d = DrawTrSurf::GetCurve(theArgv[1]);
1717   if(GC3d.IsNull())
1718   {
1719     GC2d = DrawTrSurf::GetCurve2d(theArgv[1]);
1720     if(GC2d.IsNull())
1721     {
1722       theDI << "Argument is not a 2D or 3D curve!\n";
1723       return 1;
1724     }
1725     else
1726     {
1727       theDI << theArgv[1] << " has " << aContName[GC2d->Continuity()] << " continuity.\n";
1728     }
1729   }
1730   else
1731   {
1732     theDI << theArgv[1] << " has " << aContName[GC3d->Continuity()] << " continuity.\n";
1733   }
1734
1735   return 0;
1736 }
1737
1738 //=======================================================================
1739 //function : CurveCommands
1740 //purpose  :
1741 //=======================================================================
1742 void  GeometryTest::CurveCommands(Draw_Interpretor& theCommands)
1743 {
1744
1745   static Standard_Boolean loaded = Standard_False;
1746   if (loaded) return;
1747   loaded = Standard_True;
1748
1749   DrawTrSurf::BasicCommands(theCommands);
1750
1751   const char* g;
1752
1753   g = "GEOMETRY curves creation";
1754
1755   theCommands.Add("law",
1756                   "law  name degree nbknots  knot, umult  value",
1757                   __FILE__,
1758                   polelaw,g);
1759
1760   theCommands.Add("to2d","to2d c2dname c3d [plane (XOY)]",
1761                   __FILE__,
1762                   to2d,g);
1763
1764   theCommands.Add("to3d","to3d c3dname c2d [plane (XOY)]",
1765                   __FILE__,
1766                   to3d,g);
1767
1768   theCommands.Add("gproject",
1769                   "gproject projectname curve surface [tolerance [maxdist]]\n"
1770                   "\t\t[-c continuity][-d maxdegree][-s maxsegments][-2d proj2d][-3d proj3d]\n"
1771                   "\t\t-c continuity  : set curve continuity (C0, C1, C2) for approximation\n"
1772                   "\t\t-d maxdegree   : set max possible degree of result for approximation\n"
1773                   "\t\t-s maxsegments : set max value of parametric intervals the projected curve for approximation\n"
1774                   "\t\t-2d proj2d     : set necessity of 2d results (0 or 1)\n"
1775                   "\t\t-3d proj3d     : set necessity of 3d results (0 or 1)",
1776                   __FILE__,
1777                   gproject,g);
1778
1779   theCommands.Add("project",
1780                   "project : no args to have help",
1781                   __FILE__,
1782                   project,g);
1783
1784   theCommands.Add("projonplane",
1785                   "projonplane r C3d Plane [dx dy dz] [0/1]",
1786                   __FILE__,
1787                   projonplane, g);
1788
1789   theCommands.Add("bisec",
1790                   "bisec result line/circle/point line/circle/point",
1791                   __FILE__,
1792                   bisec, g);
1793
1794   g = "GEOMETRY Curves and Surfaces modification";
1795
1796
1797   theCommands.Add("movelaw",
1798                   "movelaw name u  x  tx [ constraint = 0]",
1799                   __FILE__,
1800                   movelaw,g) ;
1801
1802
1803
1804   g = "GEOMETRY intersections";
1805
1806   theCommands.Add("intersect",
1807                   "intersect result surf1/curv1 surf2 [tolerance]\n\t\t  "
1808                   "intersect result surf1 surf2 [u1 v1 u2 v2] [U1F U1L V1F V1L U2F U2L V2F V2L] [tolerance]",
1809                   __FILE__,
1810                   intersection,g);
1811
1812   theCommands.Add("crvpoints",
1813                   "crvpoints result <curve or wire> deflection",
1814                   __FILE__,
1815                   crvpoints,g);
1816
1817   theCommands.Add("crvtpoints",
1818                   "crvtpoints result <curve or wire> deflection angular deflection - tangential deflection points",
1819                   __FILE__,
1820                   crvtpoints,g);
1821
1822   theCommands.Add("uniformAbscissa",
1823                   "uniformAbscissa Curve nbPnt",
1824                   __FILE__,
1825                   uniformAbscissa,g);
1826
1827   theCommands.Add("uniformAbscissaEl",
1828                   "uniformAbscissaEl maxR minR nbPnt",
1829                   __FILE__,  EllipsUniformAbscissa,g);
1830
1831   theCommands.Add("discrCurve",
1832     "discrCurve polyline curve params\n"
1833     "Approximates a curve by a polyline (first degree B-spline).\n"
1834     "nbPnts number - creates polylines with the number points\n"
1835     "uniform 0 | 1 - creates polyline with equal length segments",
1836     __FILE__,  discrCurve, g);
1837
1838   theCommands.Add("mypoints",
1839                   "mypoints result curv deflection",
1840                   __FILE__,
1841                   mypoints,g);
1842   theCommands.Add("surfpoints",
1843                   "surfoints result surf deflection",
1844                   __FILE__,
1845                   surfpoints,g);
1846
1847   theCommands.Add("getcurvcontinuity",
1848                   "getcurvcontinuity {curve or 2dcurve}: \n\tReturns the continuity of the given curve",
1849                   __FILE__,
1850                   GetCurveContinuity,g);
1851
1852
1853 }
1854