0024428: Implementation of LGPL license
[occt.git] / src / Approx / Approx_BSplComputeLine.gxx
1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and / or modify it
7 // under the terms of the GNU Lesser General Public version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 #include <stdio.h>
16
17 #include <Approx_ParametrizationType.hxx>
18 #include <TColStd_Array1OfReal.hxx>
19 #include <TColgp_Array1OfPnt.hxx>
20 #include <TColgp_Array1OfPnt2d.hxx>
21 #include <gp_Pnt.hxx>
22 #include <gp_Pnt2d.hxx>
23 #include <gp_Vec.hxx>
24 #include <gp_Vec2d.hxx>
25 #include <TColgp_Array1OfVec.hxx>
26 #include <TColgp_Array1OfVec2d.hxx>
27 #include <AppParCurves_Constraint.hxx>
28 #include <AppParCurves_HArray1OfConstraintCouple.hxx>
29 #include <AppParCurves_MultiPoint.hxx>
30 #include <Precision.hxx>
31 #include <math_IntegerVector.hxx>
32 #include <math_Gauss.hxx>
33 #include <math_Uzawa.hxx>
34 #include <AppParCurves_ConstraintCouple.hxx>
35 #include Approx_BSpParLeastSquareOfMyBSplGradient_hxx
36
37 #if defined( DEB ) && defined( DRAW ) && !defined( WNT )
38
39 static Standard_Boolean mydebug = Standard_False;
40
41 #include <Draw.hxx>
42 #include <Draw_Appli.hxx>
43 #include <DrawTrSurf.hxx>
44 #include <Draw_Text2D.hxx>
45 #include <Draw_Text3D.hxx>
46 #include <TColStd_Array1OfInteger.hxx>
47 #include <Geom_BSplineCurve.hxx>
48 #include <Geom2d_BSplineCurve.hxx>
49 #include <Geom_Line.hxx>
50 #include <Geom2d_Line.hxx>
51 #include <Geom_TrimmedCurve.hxx>
52 #include <Geom2d_TrimmedCurve.hxx>
53
54
55 static void DUMP(const MultiLine& Line)
56 {
57   Standard_Integer i, j, nbP2d, nbP3d, firstP, lastP;
58   firstP = LineTool::FirstPoint(Line);
59   lastP  = LineTool::LastPoint(Line);
60
61   nbP3d = LineTool::NbP3d(Line);
62   nbP2d = LineTool::NbP2d(Line);
63   Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
64   if (nbP3d == 0) mynbP3d = 1;
65   if (nbP2d == 0) mynbP2d = 1;
66   
67   TColgp_Array1OfPnt tabP(1, mynbP3d);
68   TColgp_Array1OfPnt2d tabP2d(1, mynbP2d);
69   TColgp_Array1OfVec TabV(1, mynbP3d);
70   TColgp_Array1OfVec2d TabV2d(1, mynbP2d);
71   Standard_Boolean Ok;
72   Handle(Geom_Line) L3d;
73   Handle(Geom2d_Line) L2d;
74   Handle(Geom_TrimmedCurve) L3dt;
75   Handle(Geom2d_TrimmedCurve) L2dt;
76   Handle(Draw_Text2D) T2D;
77   Handle(Draw_Text3D) T3D;
78
79   char solname[100];
80   char mytext[10];
81
82   for (i = firstP; i <= lastP; i++) {
83     if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, i, tabP, tabP2d);
84     else if (nbP2d != 0)          LineTool::Value(Line, i, tabP2d);
85     else if (nbP3d != 0)          LineTool::Value(Line, i, tabP);
86     
87     for (j = 1; j <= nbP3d; j++) {
88       sprintf(solname, "%s%i%s_%i", "p", j, "3d", i);
89       char* Temp = solname ;
90       DrawTrSurf::Set(Temp, tabP(j));
91 //      DrawTrSurf::Set(solname, tabP(j));
92       if (i == firstP || i == lastP) {
93         sprintf(mytext, "%s%i", " ", i);
94         T3D = new Draw_Text3D(tabP(j), mytext, Draw_vert);
95         dout << T3D;
96       }
97     }
98     for (j = 1; j <= nbP2d; j++) {
99       sprintf(solname, "%s%i%s_%i", "p", j, "2d", i);
100       char* Temp = solname ;
101       DrawTrSurf::Set(Temp, tabP2d(j));
102 //      DrawTrSurf::Set(solname, tabP2d(j));
103       if (i == firstP || i == lastP) {
104         sprintf(mytext, "%s%i", " ", i);
105         T2D = new Draw_Text2D(tabP2d(j), mytext, Draw_vert);
106         dout << T2D;
107       }
108     }
109
110     // le cas des tangentes aux extremites:
111     if (i == firstP || i == lastP) {
112       if (nbP3d != 0 && nbP2d != 0)
113         Ok = LineTool::Tangency(Line, i, TabV, TabV2d);
114       else if (nbP2d != 0)
115         Ok = LineTool::Tangency(Line, i, TabV2d);
116       else if (nbP3d != 0)
117         Ok = LineTool::Tangency(Line, i, TabV);
118       
119       if (Ok) {
120         for (j = 1; j <= nbP3d; j++) {
121           sprintf(solname, "%s%i%s_%i", "t", j, "3d", i);
122           L3d = new Geom_Line (tabP(j), gp_Dir(TabV(j)));
123           L3dt = new Geom_TrimmedCurve(L3d, 0.0, 0.3);
124           char* Temp = solname ;
125           DrawTrSurf::Set(Temp, L3dt);
126 //        DrawTrSurf::Set(solname, L3dt);
127         }
128         for (j = 1; j <= nbP2d; j++) {
129           sprintf(solname, "%s%i%s_%i", "t", j, "2d", i);
130           L2d = new Geom2d_Line (tabP2d(j), gp_Dir2d(TabV2d(j)));
131           L2dt = new Geom2d_TrimmedCurve(L2d, 0.0, 0.3);
132           char* Temp = solname ;
133           DrawTrSurf::Set(Temp, L2dt);
134 //        DrawTrSurf::Set(solname, L2dt);
135         }
136       }
137     }
138   }
139   dout.Flush();
140 }
141
142
143 static void DUMP(const AppParCurves_MultiBSpCurve& C)
144 {
145   static Standard_Integer nbappel = 0;
146   Standard_Integer i, j, nbP2d, nbP3d;
147   Standard_Integer nbpoles = C.NbPoles();
148   Standard_Integer deg = C.Degree();
149   const TColStd_Array1OfReal& Knots = C.Knots();
150   const TColStd_Array1OfInteger& Mults = C.Multiplicities();
151
152   Handle(Geom_BSplineCurve) BSp;
153   Handle(Geom2d_BSplineCurve) BSp2d;
154
155   TColgp_Array1OfPnt tabPoles(1, nbpoles);
156   TColgp_Array1OfPnt2d tabPoles2d(1, nbpoles);
157   char solname[100];
158
159   nbappel++;
160   for (i = 1; i <= C.NbCurves(); i++) {
161     if (C.Dimension(i) == 3) {
162       C.Curve(i, tabPoles);
163       BSp = new Geom_BSplineCurve(tabPoles, Knots, Mults, deg);
164       sprintf(solname, "%s%i%s_%i", "c", i, "3d", nbappel);
165       char* Temp = solname ;
166       DrawTrSurf::Set(Temp, BSp);
167 //      DrawTrSurf::Set(solname, BSp);
168     }
169     else {
170       C.Curve(i, tabPoles2d);
171       BSp2d = new Geom2d_BSplineCurve(tabPoles2d, Knots, Mults, deg);
172       sprintf(solname, "%s%i%s_%i", "c", i, "2d", nbappel);
173       char* Temp = solname ;
174       DrawTrSurf::Set(Temp, BSp2d);
175 //      DrawTrSurf::Set(solname, BSp2d);
176     }
177   }
178   dout.Flush();
179 }
180
181
182 #endif
183
184
185
186
187 //=======================================================================
188 //function : FirstTangencyVector
189 //purpose  : 
190 //=======================================================================
191 void Approx_BSplComputeLine::FirstTangencyVector(const MultiLine&       Line,
192                                                  const Standard_Integer index,
193                                                  math_Vector&           V)
194 const {
195
196   Standard_Integer i, j, nbP2d, nbP3d;
197   nbP3d = LineTool::NbP3d(Line);
198   nbP2d = LineTool::NbP2d(Line);
199   Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
200   if (nbP3d == 0) mynbP3d = 1;
201   if (nbP2d == 0) mynbP2d = 1;
202   Standard_Boolean Ok=Standard_False;
203   TColgp_Array1OfVec TabV(1, mynbP3d);
204   TColgp_Array1OfVec2d TabV2d(1, mynbP2d);
205  
206   if (nbP3d != 0 && nbP2d != 0)
207     Ok = LineTool::Tangency(Line, index, TabV, TabV2d);
208   else if (nbP2d != 0)
209     Ok = LineTool::Tangency(Line, index, TabV2d);
210   else if (nbP3d != 0)
211     Ok = LineTool::Tangency(Line, index, TabV);
212
213   if (Ok) {
214     if (nbP3d != 0) {
215       j = 1;
216       for (i = TabV.Lower(); i <= TabV.Upper(); i++) {
217         V(j)   = TabV(i).X();
218         V(j+1) = TabV(i).Y();
219         V(j+2) = TabV(i).Z();
220         j += 3;
221       }
222     }
223     if (nbP2d != 0) {
224       j = nbP3d*3+1;
225       for (i = TabV2d.Lower(); i <= TabV2d.Upper(); i++) {
226         V(j)   = TabV2d(i).X();
227         V(j+1) = TabV2d(i).Y();
228         j += 2;
229       }
230     }
231   }
232   else {
233
234     // recherche d un vecteur tangent par construction d une parabole:
235     AppParCurves_Constraint firstC, lastC;
236     firstC = lastC = AppParCurves_PassPoint;
237     Standard_Integer nbpoles = 3;
238     math_Vector mypar(index, index+2);
239     Parameters(Line, index, index+2, mypar);
240     Approx_BSpParLeastSquareOfMyBSplGradient 
241       LSQ(Line, index, index+2, firstC, lastC, mypar, nbpoles);
242     AppParCurves_MultiCurve C = LSQ.BezierValue();
243     
244     gp_Pnt myP;
245     gp_Vec myV;
246     gp_Pnt2d myP2d;
247     gp_Vec2d myV2d;
248     j = 1;
249     for (i = 1; i <= nbP3d; i++) {
250       C.D1(i, 0.0, myP, myV);
251       V(j)   = myV.X();
252       V(j+1) = myV.Y();
253       V(j+2) = myV.Z();
254       j += 3;
255     }
256     j = nbP3d*3+1;
257     for (i = nbP3d+1; i <= nbP3d+nbP2d; i++) {
258       C.D1(i, 0.0, myP2d, myV2d);
259       V(j)   = myV2d.X();
260       V(j+1) = myV2d.Y();
261       j += 2;
262     }
263   }
264
265 }
266
267
268 //=======================================================================
269 //function : LastTangencyVector
270 //purpose  : 
271 //=======================================================================
272 void Approx_BSplComputeLine::LastTangencyVector(const MultiLine&       Line,
273                                                 const Standard_Integer index,
274                                                 math_Vector&           V)
275 const {
276   Standard_Integer i, j, nbP2d, nbP3d;
277   nbP3d = LineTool::NbP3d(Line);
278   nbP2d = LineTool::NbP2d(Line);
279   Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
280   if (nbP3d == 0) mynbP3d = 1;
281   if (nbP2d == 0) mynbP2d = 1;
282   Standard_Boolean Ok=Standard_False;
283   TColgp_Array1OfVec TabV(1, mynbP3d);
284   TColgp_Array1OfVec2d TabV2d(1, mynbP2d);
285
286  
287   if (nbP3d != 0 && nbP2d != 0)
288     Ok = LineTool::Tangency(Line, index, TabV, TabV2d);
289   else if (nbP2d != 0)
290     Ok = LineTool::Tangency(Line, index, TabV2d);
291   else if (nbP3d != 0)
292     Ok = LineTool::Tangency(Line, index, TabV);
293
294   if (Ok) {
295     if (nbP3d != 0) {
296       j = 1;
297       for (i = TabV.Lower(); i <= TabV.Upper(); i++) {
298         V(j)   = TabV(i).X();
299         V(j+1) = TabV(i).Y();
300         V(j+2) = TabV(i).Z();
301         j += 3;
302       }
303     }
304     if (nbP2d != 0) {
305       j = nbP3d*3+1;
306       for (i = TabV2d.Lower(); i <= TabV2d.Upper(); i++) {
307         V(j)   = TabV2d(i).X();
308         V(j+1) = TabV2d(i).Y();
309         j += 2;
310       }
311     }
312   }
313   else {
314
315     // recherche d un vecteur tangent par construction d une parabole:
316     AppParCurves_Constraint firstC, lastC;
317     firstC = lastC = AppParCurves_PassPoint;
318     Standard_Integer nbpoles = 3;
319     math_Vector mypar(index-2, index);
320     Parameters(Line, index-2, index, mypar);
321     Approx_BSpParLeastSquareOfMyBSplGradient 
322       LSQ(Line, index-2, index, firstC, lastC, mypar, nbpoles);
323     AppParCurves_MultiCurve C = LSQ.BezierValue();
324     
325     gp_Pnt myP;
326     gp_Vec myV;
327     gp_Pnt2d myP2d;
328     gp_Vec2d myV2d;
329     j = 1;
330     for (i = 1; i <= nbP3d; i++) {
331       C.D1(i, 1.0, myP, myV);
332       V(j)   = myV.X();
333       V(j+1) = myV.Y();
334       V(j+2) = myV.Z();
335       j += 3;
336     }
337     j = nbP3d*3+1;
338     for (i = nbP3d+1; i <= nbP3d+nbP2d; i++) {
339       C.D1(i, 1.0, myP2d, myV2d);
340       V(j)   = myV2d.X();
341       V(j+1) = myV2d.Y();
342       j += 2;
343     }
344   }
345
346 }
347
348
349
350 //=======================================================================
351 //function : SearchFirstLambda
352 //purpose  : 
353 //=======================================================================
354 Standard_Real Approx_BSplComputeLine::
355   SearchFirstLambda(const MultiLine&            Line, 
356                     const math_Vector&          aPar,
357                     const TColStd_Array1OfReal& Theknots,
358                     const math_Vector&          V,
359                     const Standard_Integer      index) const {
360
361   // dq/dw = lambda* V = (p2-p1)/(u2-u1)
362   
363   Standard_Integer nbP2d, nbP3d;
364   gp_Pnt P1, P2;
365   gp_Pnt2d P12d, P22d;
366   nbP3d = LineTool::NbP3d(Line);
367   nbP2d = LineTool::NbP2d(Line);
368   Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
369   if (nbP3d == 0) mynbP3d = 1;
370   if (nbP2d == 0) mynbP2d = 1;
371   TColgp_Array1OfPnt tabP1(1, mynbP3d), tabP2(1, mynbP3d);
372   TColgp_Array1OfPnt2d tabP12d(1, mynbP2d), tabP22d(1, mynbP2d);
373  
374
375   if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, index, tabP1, tabP12d);
376   else if (nbP2d != 0)          LineTool::Value(Line, index, tabP12d);
377   else if (nbP3d != 0)          LineTool::Value(Line, index, tabP1);
378
379   if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, index+1, tabP2, tabP22d);
380   else if (nbP2d != 0)          LineTool::Value(Line, index+1, tabP22d);
381   else if (nbP3d != 0)          LineTool::Value(Line, index+1, tabP2);
382                                      
383
384   Standard_Real U1 = aPar(index), U2 = aPar(index+1);
385   Standard_Real lambda, S;                                        
386   Standard_Integer low = V.Lower();
387   Standard_Integer nbknots = Theknots.Length();
388  
389   if (nbP3d != 0) {
390     P1 = tabP1(1);
391     P2 = tabP2(1);
392     gp_Vec P1P2(P1, P2), myV;
393     myV.SetCoord(V(low), V(low+1), V(low+2));
394     lambda = (P1P2.Magnitude())/(myV.Magnitude()*(U2-U1));
395     S = (P1P2.Dot(myV)> 0.0) ? 1.0 : -1.0;
396   }
397   else {
398     P12d = tabP12d(1);
399     P22d = tabP22d(1);
400     gp_Vec2d P1P2(P12d, P22d), myV;
401     myV.SetCoord(V(low), V(low+1));
402     lambda = (P1P2.Magnitude())/(myV.Magnitude()*(U2-U1));
403     S = (P1P2.Dot(myV)> 0.0) ? 1.0 : -1.0;
404   }
405   return ((S*lambda)*(Theknots(2)-Theknots(1))/(Theknots(nbknots)-Theknots(1)));
406
407 }
408
409
410 //=======================================================================
411 //function : SearchLastLambda
412 //purpose  : 
413 //=======================================================================
414 Standard_Real Approx_BSplComputeLine::
415   SearchLastLambda(const MultiLine&            Line, 
416                    const math_Vector&          aPar,
417                    const TColStd_Array1OfReal& Theknots,
418                    const math_Vector&          V,
419                    const Standard_Integer      index) const
420
421 {
422   // dq/dw = lambda* V = (p2-p1)/(u2-u1)
423   
424   Standard_Integer nbP2d, nbP3d;
425   gp_Pnt P1, P2;
426   gp_Pnt2d P12d, P22d;
427   nbP3d = LineTool::NbP3d(Line);
428   nbP2d = LineTool::NbP2d(Line);
429   Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
430   if (nbP3d == 0) mynbP3d = 1;
431   if (nbP2d == 0) mynbP2d = 1;
432   TColgp_Array1OfPnt tabP(1, mynbP3d), tabP2(1, mynbP3d);
433   TColgp_Array1OfPnt2d tabP2d(1, mynbP2d), tabP22d(1, mynbP2d);
434  
435   if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, index-1, tabP, tabP2d);
436   else if (nbP2d != 0)          LineTool::Value(Line, index-1, tabP2d);
437   else if (nbP3d != 0)          LineTool::Value(Line, index-1, tabP);
438
439   if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, index, tabP2, tabP22d);
440   else if (nbP2d != 0)          LineTool::Value(Line, index, tabP22d);
441   else if (nbP3d != 0)          LineTool::Value(Line, index, tabP2);
442
443
444   Standard_Real U1 = aPar(index-1), U2 = aPar(index);
445   Standard_Real lambda, S;
446   Standard_Integer low = V.Lower();
447   Standard_Integer nbknots = Theknots.Length();
448   if (nbP3d != 0) {
449     P1 = tabP(1);
450     P2 = tabP2(1);
451     gp_Vec P1P2(P1, P2), myV;
452     myV.SetCoord(V(low), V(low+1), V(low+2));
453     lambda = (P1P2.Magnitude())/(myV.Magnitude()*(U2-U1));
454     S = (P1P2.Dot(myV)> 0.0) ? 1.0 : -1.0;
455   }
456   else {
457     P12d = tabP2d(1);
458     P22d = tabP22d(1);
459     gp_Vec2d P1P2(P12d, P22d), myV;
460     myV.SetCoord(V(low), V(low+1));
461     lambda = (P1P2.Magnitude())/(myV.Magnitude()*(U2-U1));
462     S = (P1P2.Dot(myV)> 0.0) ? 1.0 : -1.0;
463   }
464
465   return ((S*lambda)*(Theknots(nbknots)-Theknots(nbknots-1))
466          /(Theknots(nbknots)-Theknots(1)));
467 }
468
469
470
471 //=======================================================================
472 //function : Approx_BSplComputeLine
473 //purpose  : 
474 //=======================================================================
475 Approx_BSplComputeLine::Approx_BSplComputeLine
476                     (const MultiLine& Line,
477                      const math_Vector& Parameters,
478                      const Standard_Integer degreemin,
479                      const Standard_Integer degreemax,
480                      const Standard_Real Tolerance3d,
481                      const Standard_Real Tolerance2d,
482                      const Standard_Integer NbIterations,
483                      const Standard_Boolean cutting,
484                      const Standard_Boolean Squares)
485 {
486   myfirstParam = new TColStd_HArray1OfReal(Parameters.Lower(), 
487                                            Parameters.Upper());
488   for (Standard_Integer i = Parameters.Lower(); i <= Parameters.Upper(); i++) {
489     myfirstParam->SetValue(i, Parameters(i));
490   }
491   myConstraints = new AppParCurves_HArray1OfConstraintCouple(1, 2);
492   Par = Approx_IsoParametric;
493   mydegremin = degreemin;
494   mydegremax = degreemax;
495   mytol3d = Tolerance3d;
496   mytol2d = Tolerance2d;
497   mysquares = Squares;
498   mycut = cutting;
499   myitermax = NbIterations;
500   alldone = Standard_False;
501   mycont = -1;
502   myfirstC = AppParCurves_TangencyPoint;
503   mylastC = AppParCurves_TangencyPoint;
504   myhasknots = Standard_False;
505   myhasmults = Standard_False;
506   currenttol3d = currenttol2d = RealLast();
507   tolreached = Standard_False;
508   Perform(Line);
509 }
510
511
512 //=======================================================================
513 //function : Approx_BSplComputeLine
514 //purpose  : 
515 //=======================================================================
516 Approx_BSplComputeLine::Approx_BSplComputeLine
517                     (const math_Vector& Parameters,
518                      const Standard_Integer degreemin,
519                      const Standard_Integer degreemax,
520                      const Standard_Real Tolerance3d,
521                      const Standard_Real Tolerance2d,
522                      const Standard_Integer NbIterations,
523                      const Standard_Boolean cutting,
524                      const Standard_Boolean Squares)
525 {
526   myfirstParam = new TColStd_HArray1OfReal(Parameters.Lower(), 
527                                            Parameters.Upper());
528   for (Standard_Integer i = Parameters.Lower(); i <= Parameters.Upper(); i++) {
529     myfirstParam->SetValue(i, Parameters(i));
530   }
531   myfirstC = AppParCurves_TangencyPoint;
532   mylastC = AppParCurves_TangencyPoint;
533   myConstraints = new AppParCurves_HArray1OfConstraintCouple(1, 2);
534   Par = Approx_IsoParametric;
535   mydegremin = degreemin;
536   mydegremax = degreemax;
537   mytol3d = Tolerance3d;
538   mytol2d = Tolerance2d;
539   mysquares = Squares;
540   mycut = cutting;
541   myitermax = NbIterations;
542   alldone = Standard_False;
543   myhasknots = Standard_False;
544   myhasmults = Standard_False;
545   mycont = -1;
546   currenttol3d = currenttol2d = RealLast();
547   tolreached = Standard_False;
548 }
549
550 //=======================================================================
551 //function : Approx_BSplComputeLine
552 //purpose  : 
553 //=======================================================================
554 Approx_BSplComputeLine::Approx_BSplComputeLine
555                     (const Standard_Integer degreemin,
556                      const Standard_Integer degreemax,
557                      const Standard_Real Tolerance3d,
558                      const Standard_Real Tolerance2d,
559                      const Standard_Integer NbIterations,
560                      const Standard_Boolean cutting,
561                      const Approx_ParametrizationType parametrization,
562                      const Standard_Boolean Squares)
563 {
564   myConstraints = new AppParCurves_HArray1OfConstraintCouple(1, 2);
565   Par = parametrization;
566   mydegremin = degreemin;
567   mydegremax = degreemax;
568   mytol3d = Tolerance3d;
569   mytol2d = Tolerance2d;
570   mysquares = Squares;
571   mycut = cutting;
572   myitermax = NbIterations;
573   myfirstC = AppParCurves_TangencyPoint;
574   mylastC = AppParCurves_TangencyPoint;
575   alldone = Standard_False;
576   myhasknots = Standard_False;
577   myhasmults = Standard_False;
578   mycont = -1;
579   currenttol3d = currenttol2d = RealLast();
580   tolreached = Standard_False;
581 }
582
583
584 //=======================================================================
585 //function : Approx_BSplComputeLine
586 //purpose  : 
587 //=======================================================================
588 Approx_BSplComputeLine::Approx_BSplComputeLine
589                     (const MultiLine& Line,
590                      const Standard_Integer degreemin,
591                      const Standard_Integer degreemax,
592                      const Standard_Real Tolerance3d,
593                      const Standard_Real Tolerance2d,
594                      const Standard_Integer NbIterations,
595                      const Standard_Boolean cutting,
596                      const Approx_ParametrizationType parametrization,
597                      const Standard_Boolean Squares)
598 {
599   myConstraints = new AppParCurves_HArray1OfConstraintCouple(1, 2);
600   alldone = Standard_False;
601   mydegremin = degreemin;
602   mydegremax = degreemax;
603   mytol3d = Tolerance3d;
604   mytol2d = Tolerance2d;
605   mysquares = Squares;
606   mycut = cutting;
607   myitermax = NbIterations;
608   Par = parametrization;
609   myfirstC = AppParCurves_TangencyPoint;
610   mylastC = AppParCurves_TangencyPoint;
611   myhasknots = Standard_False;
612   myhasmults = Standard_False;
613   mycont = -1;
614   currenttol3d = currenttol2d = RealLast();
615   tolreached = Standard_False;
616   Perform(Line);
617 }
618
619
620
621 //=======================================================================
622 //function : Perform
623 //purpose  : 
624 //=======================================================================
625 void Approx_BSplComputeLine::Perform(const MultiLine& Line)
626 {
627
628 #if defined( DEB ) && defined( DRAW ) && !defined( WNT )
629   if (mydebug) DUMP(Line);
630 #endif
631     
632   Standard_Integer i, Thefirstpt, Thelastpt;
633   Standard_Boolean Finish = Standard_False, begin = Standard_True;
634   
635   // recherche des vraies contraintes donnees par la Line:
636   FindRealConstraints(Line);
637
638   Thefirstpt = LineTool::FirstPoint(Line);
639   Thelastpt  = LineTool::LastPoint(Line);
640   Standard_Integer myfirstpt = Thefirstpt; 
641   Standard_Integer mylastpt = Thelastpt;
642
643   AppParCurves_ConstraintCouple myCouple1(myfirstpt, realfirstC);
644   AppParCurves_ConstraintCouple myCouple2(mylastpt, reallastC);
645   myConstraints->SetValue(1, myCouple1);
646   myConstraints->SetValue(2, myCouple2);
647
648   math_Vector TheParam(Thefirstpt, Thelastpt, 0.0);
649   if(myfirstParam.IsNull()) {
650     Parameters(Line, Thefirstpt, Thelastpt, TheParam);
651   }
652   else {
653     for (i = myfirstParam->Lower(); i <= myfirstParam->Upper(); i++) {
654       TheParam(i+Thefirstpt-1) = myfirstParam->Value(i);
655     }
656   }
657
658   myParameters = new TColStd_HArray1OfReal(TheParam.Lower(), TheParam.Upper());
659   for ( i = TheParam.Lower(); i <= TheParam.Upper(); i++) {
660     myParameters->SetValue(i, TheParam(i));
661   }
662   Standard_Integer nbknots = 2;
663   Standard_Real l;
664   alldone = Standard_False;
665
666   if (!mycut) {
667
668     // cas ou on ne desire pas de noeuds supplementaires.
669     // ==================================================
670
671     if (!myhasknots) {
672       TColStd_Array1OfReal theknots(1, 2);
673       TColStd_Array1OfInteger themults(1, 2);
674       theknots(1) = 0.0;
675       theknots(2) = 1.0;
676       alldone = Compute(Line, myfirstpt, mylastpt, TheParam, theknots, themults);
677     }
678     else {
679       if (!myhasmults) {
680         TColStd_Array1OfInteger themults(1, myknots->Length());
681         alldone = Compute(Line, myfirstpt, mylastpt, TheParam, 
682                           myknots->Array1(), themults);
683       }
684       else {
685         alldone = Compute(Line, myfirstpt, mylastpt, TheParam, 
686                           myknots->Array1(), mymults->ChangeArray1());
687       }
688     }
689   }
690   else {
691
692     // cas ou on va iterer a partir de noeuds donnes par l''utilisateur 
693     // ou a partir d''une bezier.
694     // ================================================================
695
696     while (!Finish) {
697
698       currenttol3d = currenttol2d = RealLast();
699
700       if (myhasknots && begin) {
701
702         if (!myhasmults) {
703
704           // 1er cas: l''utilisateur donne des noeuds de depart mais
705           // a nous de fixer  les multiplicites  en  fonction  de  la 
706           // continuite desiree.
707           // ========================================================
708
709           TColStd_Array1OfInteger TheMults(1, myknots->Length());
710           alldone = Compute(Line, myfirstpt, mylastpt, TheParam, 
711                        myknots->Array1(), TheMults);
712         }
713         else {
714
715           // 2eme cas: l''utilisateur donne des noeuds de depart 
716           // avec leurs multiplicites.
717           // ===================================================
718
719           alldone = Compute(Line, myfirstpt, mylastpt, TheParam, 
720                        myknots->Array1(), mymults->ChangeArray1());
721         }
722         begin = Standard_False;
723       }
724
725       else {
726         
727         // 3eme cas: l''utilisateur ne donne aucun noeuds de depart
728         // ========================================================
729
730         TColStd_Array1OfReal Theknots(1, nbknots);
731         TColStd_Array1OfInteger TheMults(1, nbknots);
732         Theknots(1) = 0.0;
733         Theknots(nbknots) = 1.0;
734         for (i = 2; i <= nbknots-1; i++) {
735           
736           l = (mylastpt-myfirstpt)*Standard_Real(i-1)
737                               /Standard_Real(nbknots-1);
738           Standard_Integer ll = (Standard_Integer)(l);
739           Standard_Real a = l - ll;
740           Standard_Real p1 = TheParam(ll + myfirstpt);
741           Standard_Real p2 = TheParam(ll + 1 + myfirstpt);
742           Theknots(i) = (1.-a)*p1 + a*p2;
743         }
744         
745         alldone = Compute(Line, myfirstpt, mylastpt, TheParam, Theknots, TheMults);
746         
747       }
748       
749       if (!alldone) nbknots++;
750       else Finish = Standard_True;
751     }
752   }
753
754 #if defined( DEB ) && defined( DRAW ) && !defined( WNT )
755   if (mydebug) DUMP(TheMultiBSpCurve);
756 #endif
757
758 }
759
760
761
762
763 //=======================================================================
764 //function : Parameters
765 //purpose  : 
766 //=======================================================================
767 const TColStd_Array1OfReal& Approx_BSplComputeLine::Parameters() const
768 {
769   return myParameters->Array1();
770 }
771
772
773
774 //=======================================================================
775 //function : Value
776 //purpose  : 
777 //=======================================================================
778 const AppParCurves_MultiBSpCurve& Approx_BSplComputeLine::Value() const
779 {
780   return TheMultiBSpCurve; 
781 }
782
783 //=======================================================================
784 //function : ChangeValue
785 //purpose  : 
786 //=======================================================================
787 AppParCurves_MultiBSpCurve& Approx_BSplComputeLine::ChangeValue()
788 {
789   return TheMultiBSpCurve; 
790 }
791
792 //=======================================================================
793 //function : Parameters
794 //purpose  : 
795 //=======================================================================
796
797 void Approx_BSplComputeLine::Parameters(const MultiLine& Line, 
798                                const Standard_Integer firstP,
799                                const Standard_Integer lastP,
800                                math_Vector& TheParameters) const
801 {
802   Standard_Integer i, j, Nbp, nbP2d, nbP3d;
803   Standard_Real dist;
804   gp_Pnt P1, P2;
805   gp_Pnt2d P12d, P22d;
806   Nbp = lastP-firstP+1;
807
808   
809   if (Nbp == 2) {
810     TheParameters(firstP) = 0.0;
811     TheParameters(lastP)  = 1.0;
812   }
813   else if (Par == Approx_ChordLength || Par == Approx_Centripetal) {
814     nbP3d = LineTool::NbP3d(Line);
815     nbP2d = LineTool::NbP2d(Line);
816     Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
817     if (nbP3d == 0) mynbP3d = 1;
818     if (nbP2d == 0) mynbP2d = 1;
819
820     TheParameters(firstP) = 0.0;
821     dist = 0.0;
822     TColgp_Array1OfPnt tabP(1, mynbP3d);
823     TColgp_Array1OfPnt tabPP(1, mynbP3d);
824     TColgp_Array1OfPnt2d tabP2d(1, mynbP2d);
825     TColgp_Array1OfPnt2d tabPP2d(1, mynbP2d);
826
827     for (i = firstP+1; i <= lastP; i++) {
828       if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, i-1, tabP, tabP2d);
829       else if (nbP2d != 0)          LineTool::Value(Line, i-1, tabP2d);
830       else if (nbP3d != 0)          LineTool::Value(Line, i-1, tabP);
831
832       if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, i, tabPP, tabPP2d);
833       else if (nbP2d != 0)          LineTool::Value(Line, i, tabPP2d);
834       else if (nbP3d != 0)          LineTool::Value(Line, i, tabPP);
835       dist = 0;
836       for (j = 1; j <= nbP3d; j++) {
837         P1 = tabP(j);
838         P2 = tabPP(j);
839         dist += P2.Distance(P1);
840       }
841       for (j = 1; j <= nbP2d; j++) {
842         P12d = tabP2d(j);
843         P22d = tabPP2d(j);
844         dist += P22d.Distance(P12d);
845       }
846
847       dist = dist/(nbP3d+nbP2d);
848        
849       if(Par == Approx_ChordLength)
850         TheParameters(i) = TheParameters(i-1) + dist;
851       else {// Par == Approx_Centripetal
852         TheParameters(i) = TheParameters(i-1) + Sqrt(dist);
853       }
854     }
855     for (i = firstP; i <= lastP; i++) TheParameters(i) /= TheParameters(lastP);
856   }
857   else {
858     for (i = firstP; i <= lastP; i++) {
859       TheParameters(i) = (Standard_Real(i)-firstP)/
860                          (Standard_Real(lastP)-Standard_Real(firstP));
861     }
862   }
863 }  
864
865
866 //=======================================================================
867 //function : Compute
868 //purpose  : 
869 //=======================================================================
870 Standard_Boolean Approx_BSplComputeLine::Compute(const MultiLine& Line,
871                                              const Standard_Integer fpt,
872                                              const Standard_Integer lpt,
873                                              math_Vector&     Para,
874                                              const TColStd_Array1OfReal& Knots,
875                                              TColStd_Array1OfInteger& Mults)
876 {
877   Standard_Integer i, deg, nbpoles, multinter;
878   Standard_Boolean mydone;
879   Standard_Real Fv, TheTol3d, TheTol2d, l1, l2;
880   Standard_Integer nbp = lpt-fpt+1;
881   mylambda1 = 0.0;
882   mylambda2 = 0.0;
883
884   math_Vector aParams(Para.Lower(), Para.Upper());
885
886   for (deg = mydegremin; deg <= mydegremax; deg++) {
887
888     aParams = Para;
889
890     if (!myhasmults) {
891       Mults(Mults.Lower()) = deg+1;
892       Mults(Mults.Upper()) = deg+1;
893       nbpoles = deg+1;
894       if (mycont == -1) multinter = 1;
895       else multinter = Max(1, deg-mycont);
896       for (i = Mults.Lower()+1; i <= Mults.Upper()-1; i++) {
897         Mults(i) = multinter;
898         nbpoles += multinter;
899       }
900     }
901     else {
902       nbpoles = -deg-1;
903       for (i = Mults.Lower(); i <= Mults.Upper(); i++) {
904         nbpoles += Mults.Value(i);
905       }
906     }
907
908     Standard_Integer nbpolestocompare = nbpoles;
909     if (realfirstC == AppParCurves_TangencyPoint) nbpolestocompare++;
910     if (reallastC == AppParCurves_CurvaturePoint) nbpolestocompare++;
911     if (realfirstC == AppParCurves_TangencyPoint) nbpolestocompare++;
912     if (reallastC == AppParCurves_CurvaturePoint) nbpolestocompare++;
913     if (nbpolestocompare > nbp) {
914       Interpol(Line);
915       tolreached = Standard_True;
916       return Standard_True;
917     }
918
919     AppParCurves_MultiBSpCurve mySCU(nbpoles);
920
921     if (mysquares) {
922       Approx_BSpParLeastSquareOfMyBSplGradient SQ(Line,Knots,Mults,fpt,lpt, 
923                                            realfirstC, reallastC, aParams, nbpoles);
924       mydone = SQ.IsDone();
925       if(mydone) {
926         mySCU = SQ.BSplineValue();
927         SQ.Error(Fv, TheTol3d, TheTol2d);
928       }
929       else continue;
930     }
931     else {
932       if (nbpoles != deg+1) {
933
934         if (deg == mydegremin && (realfirstC >= AppParCurves_TangencyPoint ||
935                                   reallastC >= AppParCurves_TangencyPoint)) {
936           Approx_BSpParLeastSquareOfMyBSplGradient 
937             thefitt(Line,Knots,Mults, fpt,lpt,realfirstC,reallastC, aParams, nbpoles);
938           mylambda1 = thefitt.FirstLambda()*deg;
939           mylambda2 = thefitt.LastLambda()*deg;
940          
941         }
942         l1 = mylambda1/deg;
943         l2 = mylambda2/deg;
944
945         Approx_MyBSplGradient GRAD(Line, fpt, lpt, myConstraints, 
946                                    aParams, Knots, Mults, deg, mytol3d, 
947                                    mytol2d, myitermax, l1, l2);
948
949         mydone = GRAD.IsDone();
950         if(mydone) {
951           mySCU = GRAD.Value();
952           TheTol3d = GRAD.MaxError3d();
953           TheTol2d = GRAD.MaxError2d();
954         }
955         else continue;
956       }
957       else {
958         Approx_MyGradientbis GRAD2(Line, fpt, lpt, myConstraints, 
959                                    aParams, deg, mytol3d, 
960                                    mytol2d, myitermax);
961         mydone = GRAD2.IsDone();
962         if(mydone) {
963           if (GRAD2.Value().NbCurves() == 0)
964             continue;
965           mySCU = AppParCurves_MultiBSpCurve (GRAD2.Value(), Knots, Mults);
966           TheTol3d = GRAD2.MaxError3d();
967           TheTol2d = GRAD2.MaxError2d();
968         }
969         else continue;
970       }
971     }  
972     Standard_Boolean save = Standard_True;
973     
974     for (i = aParams.Lower(); i <= aParams.Upper(); i++) {
975       if (aParams(i) <= -0.000001 || aParams(i) >= 1.000001) {
976         save = Standard_False;
977         break;
978       }
979     }
980     
981     if (mydone) {
982       if (TheTol3d <= mytol3d && TheTol2d <= mytol2d) {
983         // Stockage de la multicurve approximee.
984         tolreached = Standard_True;
985         TheMultiBSpCurve = mySCU;
986         currenttol3d = TheTol3d;
987         currenttol2d = TheTol2d;
988         if (save) {
989           for (i = aParams.Lower(); i <= aParams.Upper(); i++) {
990             myParameters->SetValue(i, aParams(i));
991           }
992         }
993         return Standard_True;
994       }
995     }
996
997     if (TheTol3d <= currenttol3d && TheTol2d <= currenttol2d) {  
998       TheMultiBSpCurve = mySCU;
999       currenttol3d = TheTol3d;
1000       currenttol2d = TheTol2d;
1001       if (save) {
1002         for (i = aParams.Lower(); i <= aParams.Upper(); i++) {
1003           myParameters->SetValue(i, aParams(i));
1004         }
1005       }
1006     }
1007
1008   }
1009
1010   return Standard_False;
1011 }
1012
1013
1014
1015 //=======================================================================
1016 //function : SetParameters
1017 //purpose  : 
1018 //=======================================================================
1019 void Approx_BSplComputeLine::SetParameters(const math_Vector& ThePar)
1020
1021   myfirstParam = new TColStd_HArray1OfReal(ThePar.Lower(), 
1022                                            ThePar.Upper());
1023   for (Standard_Integer i = ThePar.Lower(); i <= ThePar.Upper(); i++) {
1024     myfirstParam->SetValue(i, ThePar(i));
1025   }
1026 }
1027
1028
1029 //=======================================================================
1030 //function : SetKnots
1031 //purpose  : 
1032 //=======================================================================
1033 void Approx_BSplComputeLine::SetKnots(const TColStd_Array1OfReal& Knots)
1034 {
1035   myhasknots = Standard_True;
1036   myknots = new TColStd_HArray1OfReal(Knots.Lower(), Knots.Upper());
1037   for (Standard_Integer i = Knots.Lower(); i <= Knots.Upper(); i++) {
1038     myknots->SetValue(i, Knots(i));
1039   }
1040 }
1041
1042
1043 //=======================================================================
1044 //function : SetKnotsAndMultiplicities
1045 //purpose  : 
1046 //=======================================================================
1047 void Approx_BSplComputeLine::SetKnotsAndMultiplicities
1048   (const TColStd_Array1OfReal&    Knots,
1049    const TColStd_Array1OfInteger& Mults)
1050 {
1051   myhasknots = Standard_True;
1052   myhasmults = Standard_True;
1053   Standard_Integer i;
1054   myknots = new TColStd_HArray1OfReal(Knots.Lower(), Knots.Upper());
1055   for (i = Knots.Lower(); i <= Knots.Upper(); i++) {
1056     myknots->SetValue(i, Knots(i));
1057   }
1058   mymults = new TColStd_HArray1OfInteger(Mults.Lower(), Mults.Upper());
1059   for (i = Mults.Lower(); i <= Mults.Upper(); i++) {
1060     mymults->SetValue(i, Mults(i));
1061   }
1062 }
1063
1064 //=======================================================================
1065 //function : Init
1066 //purpose  : 
1067 //=======================================================================
1068 void Approx_BSplComputeLine::Init(const Standard_Integer degreemin,
1069                               const Standard_Integer degreemax,
1070                               const Standard_Real Tolerance3d,
1071                               const Standard_Real Tolerance2d,
1072                               const Standard_Integer NbIterations,
1073                               const Standard_Boolean cutting,
1074                               const Approx_ParametrizationType parametrization,
1075                               const Standard_Boolean Squares)
1076 {
1077   mydegremin = degreemin;
1078   mydegremax = degreemax;
1079   mytol3d = Tolerance3d;
1080   mytol2d = Tolerance2d;
1081   Par = parametrization;
1082   mysquares = Squares;
1083   mycut = cutting;
1084   myitermax = NbIterations;
1085 }
1086
1087
1088
1089 //=======================================================================
1090 //function : SetDegrees
1091 //purpose  : 
1092 //=======================================================================
1093 void Approx_BSplComputeLine::SetDegrees(const Standard_Integer degreemin,
1094                                     const Standard_Integer degreemax)
1095 {
1096   mydegremin = degreemin;
1097   mydegremax = degreemax;
1098 }
1099
1100
1101 //=======================================================================
1102 //function : SetTolerances
1103 //purpose  : 
1104 //=======================================================================
1105 void Approx_BSplComputeLine::SetTolerances(const Standard_Real Tolerance3d,
1106                                        const Standard_Real Tolerance2d)
1107 {
1108   mytol3d = Tolerance3d;
1109   mytol2d = Tolerance2d;
1110 }
1111
1112
1113 //=======================================================================
1114 //function : SetConstraints
1115 //purpose  : 
1116 //=======================================================================
1117 void Approx_BSplComputeLine::SetConstraints(const AppParCurves_Constraint FirstC,
1118                                         const AppParCurves_Constraint LastC)
1119 {
1120   myfirstC = FirstC;
1121   mylastC = LastC;
1122 }
1123
1124
1125
1126 //=======================================================================
1127 //function : IsAllApproximated
1128 //purpose  : 
1129 //=======================================================================
1130 Standard_Boolean Approx_BSplComputeLine::IsAllApproximated() const 
1131 {
1132   return alldone;
1133 }
1134
1135 //=======================================================================
1136 //function : IsToleranceReached
1137 //purpose  : 
1138 //=======================================================================
1139 Standard_Boolean Approx_BSplComputeLine::IsToleranceReached() const 
1140 {
1141   return tolreached;
1142 }
1143
1144 //=======================================================================
1145 //function : Error
1146 //purpose  : 
1147 //=======================================================================
1148 void Approx_BSplComputeLine::Error(Standard_Real& tol3d,
1149                                    Standard_Real& tol2d) const
1150 {
1151   tol3d = currenttol3d;
1152   tol2d = currenttol2d;
1153 }
1154
1155
1156
1157 //=======================================================================
1158 //function : SetContinuity
1159 //purpose  : 
1160 //=======================================================================
1161 void Approx_BSplComputeLine::SetContinuity(const Standard_Integer C)
1162 {
1163   mycont = C;
1164 }
1165
1166
1167
1168 //=======================================================================
1169 //function : FindRealConstraints
1170 //purpose  : 
1171 //=======================================================================
1172 void Approx_BSplComputeLine::FindRealConstraints(const MultiLine& Line)
1173 {
1174   realfirstC = myfirstC;
1175   reallastC = mylastC;
1176   Standard_Integer nbP2d, nbP3d;
1177   nbP3d = LineTool::NbP3d(Line);
1178   nbP2d = LineTool::NbP2d(Line);
1179   Standard_Boolean Ok=Standard_False;
1180   TColgp_Array1OfVec TabV(1, Max(1, nbP3d));
1181   TColgp_Array1OfVec2d TabV2d(1, Max(1, nbP2d));
1182   Standard_Integer Thefirstpt = LineTool::FirstPoint(Line);
1183   Standard_Integer Thelastpt  = LineTool::LastPoint(Line);
1184
1185   if (myfirstC >= AppParCurves_TangencyPoint) {
1186     
1187     if (nbP3d != 0 && nbP2d != 0)
1188       Ok = LineTool::Tangency(Line, Thefirstpt, TabV, TabV2d);
1189     else if (nbP2d != 0)
1190       Ok = LineTool::Tangency(Line, Thefirstpt, TabV2d);
1191     else if (nbP3d != 0)
1192       Ok = LineTool::Tangency(Line, Thefirstpt, TabV);
1193
1194     realfirstC = AppParCurves_PassPoint;
1195     if (Ok) {
1196       realfirstC = AppParCurves_TangencyPoint;
1197       if (myfirstC == AppParCurves_CurvaturePoint) {
1198         if (nbP3d != 0 && nbP2d != 0)
1199           Ok = LineTool::Tangency(Line, Thefirstpt, TabV, TabV2d);
1200         else if (nbP2d != 0)
1201           Ok = LineTool::Tangency(Line, Thefirstpt, TabV2d);
1202         else if (nbP3d != 0)
1203           Ok = LineTool::Tangency(Line, Thefirstpt, TabV);
1204         if (Ok) {
1205           realfirstC = AppParCurves_CurvaturePoint;
1206         }
1207       }
1208     }
1209   }
1210
1211
1212   if (mylastC >= AppParCurves_TangencyPoint) {
1213     
1214     if (nbP3d != 0 && nbP2d != 0)
1215       Ok = LineTool::Tangency(Line, Thelastpt, TabV, TabV2d);
1216     else if (nbP2d != 0)
1217       Ok = LineTool::Tangency(Line, Thelastpt, TabV2d);
1218     else if (nbP3d != 0)
1219       Ok = LineTool::Tangency(Line, Thelastpt, TabV);
1220
1221     reallastC = AppParCurves_PassPoint;
1222     if (Ok) {
1223       reallastC = AppParCurves_TangencyPoint;
1224       if (mylastC == AppParCurves_CurvaturePoint) {
1225         if (nbP3d != 0 && nbP2d != 0)
1226           Ok = LineTool::Tangency(Line, Thelastpt, TabV, TabV2d);
1227         else if (nbP2d != 0)
1228           Ok = LineTool::Tangency(Line, Thelastpt, TabV2d);
1229         else if (nbP3d != 0)
1230           Ok = LineTool::Tangency(Line, Thelastpt, TabV);
1231         if (Ok) {
1232           reallastC = AppParCurves_CurvaturePoint;
1233         }
1234       }
1235     }
1236   }
1237
1238 }
1239
1240
1241
1242
1243 //=======================================================================
1244 //function : Interpol
1245 //purpose  : 
1246 //=======================================================================
1247 void Approx_BSplComputeLine::Interpol(const MultiLine& Line) 
1248 {
1249   Standard_Integer i, Thefirstpt, Thelastpt, deg = 3;
1250   mycont = 2;
1251   Thefirstpt = LineTool::FirstPoint(Line);
1252   Thelastpt  = LineTool::LastPoint(Line);
1253   math_Vector TheParam(Thefirstpt, Thelastpt, 0.0);
1254   //Par = Approx_ChordLength;
1255   if(myfirstParam.IsNull()) {
1256     Parameters(Line, Thefirstpt, Thelastpt, TheParam);
1257   }
1258   else {
1259     for (i = myfirstParam->Lower(); i <= myfirstParam->Upper(); i++) {
1260       TheParam(i+Thefirstpt-1) = myfirstParam->Value(i);
1261     }
1262   }
1263   AppParCurves_Constraint Cons = AppParCurves_TangencyPoint;
1264   Standard_Real lambda1, lambda2;
1265   Standard_Real Fv;
1266
1267   // Recherche du nombre de noeuds.
1268   Standard_Integer nbknots, nbpoles, nbpoints;
1269   nbpoints = Thelastpt - Thefirstpt + 1;
1270
1271   if (nbpoints == 2) {
1272     Cons = AppParCurves_NoConstraint;
1273     Standard_Integer mydeg = 1;
1274     Approx_BSpParLeastSquareOfMyBSplGradient 
1275       LSQ(Line, Thefirstpt, Thelastpt, Cons, Cons, TheParam, mydeg+1);
1276     alldone = LSQ.IsDone();
1277     TColStd_Array1OfReal TheKnots(1, 2);
1278     TColStd_Array1OfInteger TheMults(1, 2);
1279     TheKnots(1) = TheParam(Thefirstpt);  TheKnots(2) = TheParam(Thelastpt);
1280     TheMults(1) = TheMults(2) = mydeg+1;
1281     TheMultiBSpCurve = 
1282       AppParCurves_MultiBSpCurve (LSQ.BezierValue(),TheKnots,TheMults);
1283     LSQ.Error(Fv, currenttol3d, currenttol2d);
1284
1285   }
1286   else {
1287     nbpoles = nbpoints + 2;
1288     nbknots = nbpoints;
1289     
1290     // Resolution:
1291     TColStd_Array1OfReal Theknots(1, nbknots);
1292     Theknots(1) = TheParam(Thefirstpt);  
1293     Theknots(nbknots) = TheParam(Thelastpt);
1294     TColStd_Array1OfInteger TheMults(1, nbknots);
1295     TheMults(1) = deg+1;
1296     TheMults(nbknots) = deg+1;
1297     
1298     Standard_Integer low = TheParam.Lower();
1299     for (i = 2; i <= nbknots-1; i++) {
1300       Theknots(i) = TheParam(i+low-1);
1301       TheMults(i) = 1;  
1302     }
1303
1304
1305     Standard_Integer nbP = 3*LineTool::NbP3d(Line)+ 2*LineTool::NbP2d(Line);
1306     math_Vector V1(1, nbP), V2(1, nbP);
1307
1308     if (nbpoints == 3 || nbpoints == 4) {
1309       FirstTangencyVector(Line, Thefirstpt, V1);
1310       lambda1 = SearchFirstLambda(Line, TheParam, Theknots, V1, Thefirstpt);
1311       
1312       LastTangencyVector(Line, Thelastpt, V2);
1313       lambda2 = SearchLastLambda(Line, TheParam, Theknots, V2, Thelastpt);
1314     }
1315
1316     else {
1317       Standard_Integer nnpol, nnp = Min(nbpoints, 9);
1318       nnpol = nnp;
1319       Standard_Integer lastp = Min(Thelastpt, Thefirstpt + nnp-1);
1320       Standard_Real U;
1321       Approx_BSpParLeastSquareOfMyBSplGradient 
1322         SQ1(Line, Thefirstpt, lastp, Cons, Cons, nnpol);
1323       
1324       math_Vector P1(Thefirstpt, lastp);
1325       for (i=Thefirstpt; i <=lastp; i++) P1(i) = TheParam(i);
1326       SQ1.Perform(P1);
1327       const AppParCurves_MultiCurve& C1 = SQ1.BezierValue();
1328       U = 0.0;
1329       TangencyVector(Line, C1, U, V1);
1330       
1331       Standard_Integer firstp = Max(Thefirstpt, Thelastpt - nnp+1);
1332       
1333       
1334       if (firstp == Thefirstpt && lastp == Thelastpt) {
1335         U = 1.0;
1336         TangencyVector (Line, C1, U, V2);
1337       }
1338       else {
1339         Approx_BSpParLeastSquareOfMyBSplGradient 
1340           SQ2(Line, firstp, Thelastpt, Cons, Cons, nnpol);
1341         
1342         math_Vector P2(firstp, Thelastpt);
1343         for (i=firstp; i <=Thelastpt; i++) P2(i) = TheParam(i);
1344         SQ2.Perform(P2);
1345         const AppParCurves_MultiCurve& C2 = SQ2.BezierValue();
1346         U = 1.0;
1347         TangencyVector(Line, C2, U, V2);
1348       }
1349       
1350       
1351       lambda1 = 1./deg;
1352       lambda1 = lambda1*(Theknots(2)-Theknots(1))
1353         /(Theknots(nbknots)-Theknots(1));
1354       lambda2 = 1./deg;
1355       lambda2 = lambda2*(Theknots(nbknots)-Theknots(nbknots-1))
1356         /(Theknots(nbknots)-Theknots(1));
1357       
1358     }
1359
1360     Approx_BSpParLeastSquareOfMyBSplGradient 
1361       SQ(Line, Theknots,TheMults,Thefirstpt, Thelastpt, 
1362          Cons, Cons, nbpoles);
1363     
1364     lambda1 = lambda1/deg;
1365     lambda2 = lambda2/deg;
1366     SQ.Perform(TheParam, V1, V2, lambda1, lambda2);
1367     alldone = SQ.IsDone();
1368     TheMultiBSpCurve = SQ.BSplineValue();
1369     SQ.Error(Fv, currenttol3d, currenttol2d);
1370     tolreached = Standard_True;
1371   }
1372   myParameters = new TColStd_HArray1OfReal(TheParam.Lower(), TheParam.Upper());
1373   for (i = TheParam.Lower(); i <= TheParam.Upper(); i++) {
1374     myParameters->SetValue(i, TheParam(i));
1375   }
1376 }
1377
1378
1379 //=======================================================================
1380 //function : TangencyVector
1381 //purpose  : 
1382 //=======================================================================
1383 void Approx_BSplComputeLine::TangencyVector(
1384        const MultiLine&               Line,
1385        const AppParCurves_MultiCurve& C,
1386        const Standard_Real            U,
1387        math_Vector&                   V) const 
1388 {
1389
1390   Standard_Integer i, j, nbP2d, nbP3d;
1391   nbP3d = LineTool::NbP3d(Line);
1392   nbP2d = LineTool::NbP2d(Line);
1393     
1394   gp_Pnt myP;
1395   gp_Vec myV;
1396   gp_Pnt2d myP2d;
1397   gp_Vec2d myV2d;
1398   j = 1;
1399   for (i = 1; i <= nbP3d; i++) {
1400     C.D1(i, U, myP, myV);
1401     V(j)   = myV.X();
1402     V(j+1) = myV.Y();
1403     V(j+2) = myV.Z();
1404     j += 3;
1405   }
1406   j = nbP3d*3+1;
1407   for (i = nbP3d+1; i <= nbP3d+nbP2d; i++) {
1408     C.D1(i, U, myP2d, myV2d);
1409     V(j)   = myV2d.X();
1410     V(j+1) = myV2d.Y();
1411     j += 2;
1412   }
1413
1414 }
1415