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