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