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