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