0027772: Foundation Classes - define Standard_Boolean using C++ type "bool" instead...
[occt.git] / 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_Boolean 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