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