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