0022312: Translation of french commentaries in OCCT files
[occt.git] / src / Approx / Approx_SameParameter.cxx
1 // File:        Approx_SameParameter.cxx
2 // Created:     Tue Jun  6 09:51:17 1995
3 // Author:      Xavier BENVENISTE
4 //              <xab@nonox>
5
6 //  Modified by skv - Wed Jun  2 11:49:59 2004 OCC5898
7
8 #include <Approx_SameParameter.ixx>
9 #include <TColStd_Array1OfReal.hxx>
10 #include <BSplCLib.hxx>
11 #include <Adaptor3d_CurveOnSurface.hxx>
12 #include <Geom2dAdaptor_Curve.hxx>
13 #include <Geom2dAdaptor_HCurve.hxx>
14 #include <GeomAdaptor_Curve.hxx>
15 #include <GeomAdaptor_HCurve.hxx>
16 #include <GeomAdaptor_Surface.hxx>
17 #include <GeomAdaptor_HSurface.hxx>
18 //#include <GCPnts_UniformDeflection.hxx>
19 #include <GCPnts_QuasiUniformDeflection.hxx>
20 #include <Extrema_LocateExtPC.hxx>
21 #include <AdvApprox_ApproxAFunction.hxx>
22 #include <GeomLib_MakeCurvefromApprox.hxx>
23 #include <Precision.hxx>
24
25 #define MAX_ARRAY_SIZE 1000 // IFV, Jan 2000
26
27 #ifdef DEB
28 #ifdef DRAW
29 #include <DrawTrSurf.hxx>
30 #endif
31 #include <Geom2d_BSplineCurve.hxx>
32 #include <stdio.h>
33 static Standard_Boolean Voir     = Standard_False;
34 static Standard_Boolean AffichFw = Standard_False;
35 static Standard_Integer NbCurve = 0;
36 #endif
37 //
38 //   allows testing if Extrema produces correct results/
39
40
41 static void ProjectPointOnCurve(const Standard_Real      InitValue,
42                                 const gp_Pnt             APoint,
43                                 const Standard_Real      Tolerance,
44                                 const Standard_Integer   NumIteration,
45                                 const Adaptor3d_Curve&     Curve,
46                                 Standard_Boolean&        Status,
47                                 Standard_Real&           Result)
48 {
49   Standard_Integer num_iter = 0,
50   not_done = 1,
51   ii ;
52   
53   gp_Pnt a_point ;
54   gp_Vec   vector,
55   d1,
56   d2 ;
57   Standard_Real func,
58   func_derivative,
59   param = InitValue ;
60   Status = Standard_False ;
61   Standard_Real Toler = 1.0e-12;
62   do {
63     num_iter += 1 ;
64     Curve.D2(param,
65               a_point,
66              d1,
67              d2) ;
68     for (ii = 1 ; ii <= 3 ; ii++) {
69       vector.SetCoord(ii, APoint.Coord(ii) - a_point.Coord(ii)) ;
70     }
71     func = vector.Dot(d1) ;
72     func_derivative = vector.Dot(d2) ;
73     func_derivative -= d1.Dot(d1) ;
74     if ( Abs(func) < Tolerance * d1.Magnitude()) {
75       not_done = 0 ;
76       Status = Standard_True ;
77     }
78     else
79       { // fixing a bug PRO18577 : avoid divizion by zero
80         if( Abs(func_derivative) > Toler )  {
81           param -= func / func_derivative ;
82         }
83         param = Max(param,Curve.FirstParameter()) ;
84         param = Min(param,Curve.LastParameter())  ;
85         Status = Standard_True ;
86       }
87   } 
88   while (not_done && num_iter <= NumIteration) ;
89   Result = param ;
90 }  
91      
92
93
94 //=======================================================================
95 //class : Approx_SameParameter_Evaluator
96 //purpose  : 
97 //=======================================================================
98
99 class Approx_SameParameter_Evaluator : public AdvApprox_EvaluatorFunction
100 {
101  public:
102   Approx_SameParameter_Evaluator (const TColStd_Array1OfReal& theFlatKnots, 
103                                   const TColStd_Array1OfReal& thePoles,
104                                   const Handle(Adaptor2d_HCurve2d)& theHCurve2d) 
105     : FlatKnots(theFlatKnots), Poles(thePoles), HCurve2d(theHCurve2d) {}
106
107   virtual void Evaluate (Standard_Integer *Dimension,
108                          Standard_Real     StartEnd[2],
109                          Standard_Real    *Parameter,
110                          Standard_Integer *DerivativeRequest,
111                          Standard_Real    *Result, // [Dimension]
112                          Standard_Integer *ErrorCode);
113   
114  private:
115   const TColStd_Array1OfReal& FlatKnots;
116   const TColStd_Array1OfReal& Poles;
117   Handle(Adaptor2d_HCurve2d) HCurve2d;
118 };
119
120 void Approx_SameParameter_Evaluator::Evaluate (Standard_Integer *,/*Dimension*/
121                                                Standard_Real    /*StartEnd*/[2],
122                                                Standard_Real    *Parameter,
123                                                Standard_Integer *DerivativeRequest,
124                                                Standard_Real    *Result,
125                                                Standard_Integer *ReturnCode) 
126
127   gp_Pnt2d Point ;
128   gp_Vec2d Vector ;
129   Standard_Integer extrap_mode[2] ;
130   extrap_mode[0] = extrap_mode[1] = 3;
131   Standard_Real eval_result[2] ;
132   Standard_Real *PolesArray =
133     (Standard_Real *) &Poles(Poles.Lower()) ;
134   //
135   // evaluate the 1D bspline that represents the change in parameterization
136   //
137   BSplCLib::Eval(*Parameter,
138                  Standard_False,
139                  *DerivativeRequest,
140                  extrap_mode[0],
141                  3,
142                  FlatKnots,
143                  1,
144                  PolesArray[0],
145                  eval_result[0]) ;
146   
147   
148   if (*DerivativeRequest == 0){
149     HCurve2d->D0(eval_result[0],Point);
150     Point.Coord(Result[0],Result[1]);
151   }
152   else if (*DerivativeRequest == 1){
153     HCurve2d->D1(eval_result[0], Point, Vector);
154     Vector.Multiply(eval_result[1]);
155     Vector.Coord(Result[0],Result[1]);
156   }
157   ReturnCode[0] = 0 ;
158 }
159
160 static Standard_Real ComputeTolReached(const Handle(Adaptor3d_HCurve)& c3d,
161                                        const Adaptor3d_CurveOnSurface& cons,
162                                        const Standard_Integer        nbp)
163 {
164   Standard_Real d2 = 0.;
165   Standard_Integer nn = nbp;
166   Standard_Real unsurnn = 1./nn;
167   Standard_Real first = c3d->FirstParameter();
168   Standard_Real last  = c3d->LastParameter();
169   for(Standard_Integer i = 0; i <= nn; i++){
170     Standard_Real t = unsurnn*i;
171     Standard_Real u = first*(1.-t) + last*t;
172     gp_Pnt Pc3d = c3d->Value(u);
173     gp_Pnt Pcons = cons.Value(u);
174     if (Precision::IsInfinite(Pcons.X()) ||
175         Precision::IsInfinite(Pcons.Y()) ||
176         Precision::IsInfinite(Pcons.Z())) {
177       d2=Precision::Infinite();
178       break;
179     }
180     Standard_Real temp = Pc3d.SquareDistance(Pcons);
181     if(temp > d2) d2 = temp;
182   }
183   d2 = 1.5*sqrt(d2);
184   if(d2<1.e-7) d2 = 1.e-7;
185   return d2;
186 }
187
188 static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots, 
189                               const TColStd_Array1OfReal& Poles,
190                               const Standard_Integer nbp,
191                               const TColStd_Array1OfReal& pc3d,
192 //                            const TColStd_Array1OfReal& pcons,
193                               const TColStd_Array1OfReal& ,
194                               const Handle(Adaptor3d_HCurve)& c3d,
195                               const Adaptor3d_CurveOnSurface& cons,
196                               Standard_Real& tol,
197                               const Standard_Real oldtol)
198 {
199   Standard_Real d = tol;
200   Standard_Integer extrap_mode[2] ;
201   extrap_mode[0] = extrap_mode[1] = 3;
202   Standard_Integer i;
203 #ifdef DEB
204   if (Voir) {
205     cout<<endl;
206     cout<<"Control the change of variable : "<<endl;
207     cout<<"yawn mesured by projection : "<<d<<endl;
208     cout<<"Number of points : "<<nbp<<endl;
209   }
210 #endif
211 #if 0
212   Standard_Real glis = 0., dglis = 0.;
213   for(i = 1; i <= nbp; i++){
214     Standard_Real tc3d = pc3d(i);
215     gp_Pnt Pc3d = c3d->Value(tc3d);
216     Standard_Real tcons;
217     BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
218                    3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
219     gp_Pnt Pcons = cons.Value(tcons);
220     Standard_Real temp = Pc3d.SquareDistance(Pcons);
221     if(temp >= dglis) dglis = temp;
222     temp = Abs(tcons-pcons(i));
223     if(temp >= glis) glis = temp;
224   }
225   dglis = sqrt(dglis);
226 #ifdef DEB
227   if ( Voir) {
228     cout<<"shift of parameter to the imposed points : "<<glis<<endl;
229     cout<<"shift distance at the imposed points : "<<dglis<<endl;
230   }
231 #endif
232   dglis = 0.;
233   for(i = 1; i < nbp; i++){
234     Standard_Real tc3d = 0.5*(pc3d(i)+pc3d(i+1));
235     gp_Pnt Pc3d = c3d->Value(tc3d);
236     Standard_Real tcons;
237     BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
238                    3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
239     gp_Pnt Pcons = cons.Value(tcons);
240     Standard_Real temp = Pc3d.SquareDistance(Pcons);
241     if(temp >= dglis) dglis = temp;
242   }
243   dglis = sqrt(dglis);
244 #ifdef DEB
245   if (Voir)
246     cout<<"distance de glissement en milieu d intervals : "<<dglis<<endl;
247 #endif
248 #endif
249
250   Standard_Real d2 = 0.;
251   Standard_Integer nn = 2*nbp;
252   Standard_Real unsurnn = 1./nn;
253 //  Modified by skv - Wed Jun  2 11:49:59 2004 OCC5898 Begin
254 // Correction of the interval of valid values. This condition has no sensible
255 // grounds. But it is better then the old one (which is commented out) because
256 // it fixes the bug OCC5898. To develop more or less sensible criterion it is
257 // necessary to deeply investigate this problem which is not possible in frames
258 // of debugging.
259
260 //   Standard_Real firstborne= 2*pc3d(1)-pc3d(nbp);
261 //   Standard_Real lastborne= 2*pc3d(nbp)-pc3d(1);
262   Standard_Real firstborne= 3.*pc3d(1)   - 2.*pc3d(nbp);
263   Standard_Real lastborne = 3.*pc3d(nbp) - 2.*pc3d(1);
264 //  Modified by skv - Wed Jun  2 11:50:03 2004 OCC5898 End
265   for(i = 0; i <= nn; i++){
266     Standard_Real t = unsurnn*i;
267     Standard_Real tc3d = pc3d(1)*(1.-t) + pc3d(nbp)*t;
268     gp_Pnt Pc3d = c3d->Value(tc3d);
269     Standard_Real tcons;
270     BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
271                    3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
272     if (tcons < firstborne || tcons > lastborne) {
273       tol=Precision::Infinite();
274       return Standard_False;
275     }
276     gp_Pnt Pcons = cons.Value(tcons);
277     Standard_Real temp = Pc3d.SquareDistance(Pcons);
278     if(temp > d2) d2 = temp;
279   }
280   tol = sqrt(d2);
281 #ifdef DEB
282   if (Voir)
283     cout<<"distance max on "<<nn<<" points : "<<tol<<endl<<endl;
284 #endif
285   return ((tol <= d) || (tol > 0.8 * oldtol));
286 }
287
288
289 //=======================================================================
290 //function : Approx_SameParameter
291 //purpose  : 
292 //=======================================================================
293
294 Approx_SameParameter::Approx_SameParameter(const Handle(Geom_Curve)&   C3D,
295                                            const Handle(Geom2d_Curve)& C2D,
296                                            const Handle(Geom_Surface)& S,
297                                            const Standard_Real         Tol):
298  mySameParameter(Standard_True), myDone(Standard_False)
299 {
300   myHCurve2d = new Geom2dAdaptor_HCurve(C2D);
301   myC3d      = new GeomAdaptor_HCurve(C3D);
302   mySurf     = new GeomAdaptor_HSurface(S);
303   Build(Tol);
304 }
305
306
307 //=======================================================================
308 //function : Approx_SameParameter
309 //purpose  : 
310 //=======================================================================
311
312 Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)&   C3D,
313                                            const Handle(Geom2d_Curve)&     C2D,
314                                            const Handle(Adaptor3d_HSurface)& S,
315                                            const Standard_Real            Tol):
316  mySameParameter(Standard_True), myDone(Standard_False)
317 {
318   myC3d = C3D;
319   mySurf = S;
320   myHCurve2d = new Geom2dAdaptor_HCurve(C2D);
321   Build(Tol);
322 }
323
324
325 //=======================================================================
326 //function : Approx_SameParameter
327 //purpose  : 
328 //=======================================================================
329
330 Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)&   C3D,
331                                            const Handle(Adaptor2d_HCurve2d)& C2D,
332                                            const Handle(Adaptor3d_HSurface)& S,
333                                            const Standard_Real            Tol):
334  mySameParameter(Standard_True), myDone(Standard_False)
335 {
336   myC3d = C3D;
337   mySurf = S;
338   myHCurve2d = C2D;
339   Build(Tol);
340 }
341
342
343 //=======================================================================
344 //function : Build
345 //purpose  : 
346 //=======================================================================
347
348 void Approx_SameParameter::Build(const Standard_Real Tolerance)
349 {
350   Standard_Integer ii ;
351   Adaptor3d_CurveOnSurface CurveOnSurface(myHCurve2d,mySurf);
352   Standard_Real fcons = CurveOnSurface.FirstParameter();
353   Standard_Real lcons = CurveOnSurface.LastParameter();
354   Standard_Real fc3d = myC3d->FirstParameter();
355   Standard_Real lc3d = myC3d->LastParameter();
356
357   GeomAbs_Shape Continuity = myHCurve2d->Continuity();
358
359   if(Continuity > GeomAbs_C1) Continuity = GeomAbs_C1;
360
361   //Control tangents at the extremities to know if the
362   //reparametring is possible and calculate the tangents 
363   //at the extremities of the function of change of variable.
364   Standard_Real tangent[2];
365   gp_Pnt Pcons,Pc3d;
366   gp_Vec Vcons,Vc3d;
367
368   Standard_Real Tol = Tolerance;
369   Standard_Real Tol2 = Tol * Tol;
370   Standard_Real Tolp = myC3d->Resolution(Tol), deltamin = 50*Tolp;
371
372   Standard_Real besttol2 = Tol2;
373   Standard_Boolean extrok = 0;
374
375   extrok = 1;
376   CurveOnSurface.D1(fcons,Pcons,Vcons);
377   myC3d->D1(fc3d,Pc3d,Vc3d);
378   Standard_Real dist2 = Pcons.SquareDistance(Pc3d);
379   Standard_Real dmax2 = dist2;
380
381   Standard_Real magVcons = Vcons.Magnitude();
382   if (magVcons > 1.e-12){
383     tangent[0] = Vc3d.Magnitude() / magVcons;
384   }
385   else extrok = 0;
386
387   CurveOnSurface.D1(lcons,Pcons,Vcons);
388   myC3d->D1(lc3d,Pc3d,Vc3d);
389   dist2 = Pcons.SquareDistance(Pc3d);
390
391   if(dist2 > dmax2) dmax2 = dist2;
392   magVcons = Vcons.Magnitude();
393   if (magVcons > 1.e-12){
394     tangent[1] = Vc3d.Magnitude() / magVcons;
395   }
396   else extrok = 0;
397
398
399   if(dmax2 > besttol2) besttol2 = dmax2;
400
401   //Take a multiple of the sample pof CheckShape,
402   //at least the control points will be correct. No comment!!!
403
404   Standard_Integer NCONTROL = 22;
405 #ifdef DEB
406   Standard_Integer nbcoups = 0;
407 #endif
408   
409   Standard_Boolean interpolok = 0;
410   Standard_Real tolsov = 1.e200;
411   //Take parameters with  constant step on the curve on surface
412   //and on curve 3d.
413   Standard_Real deltacons = lcons - fcons;
414   deltacons /= (NCONTROL);
415   Standard_Real deltac3d = lc3d - fc3d;
416   deltac3d /= (NCONTROL);
417
418   Standard_Real wcons = fcons;
419   Standard_Real wc3d  = fc3d;
420   
421   Standard_Real qpcons[MAX_ARRAY_SIZE], qnewpcons[MAX_ARRAY_SIZE], 
422                 qpc3d[MAX_ARRAY_SIZE], qnewpc3d[MAX_ARRAY_SIZE];
423   Standard_Real * pcons = qpcons; Standard_Real * newpcons = qnewpcons;
424   Standard_Real * pc3d = qpc3d; Standard_Real * newpc3d = qnewpc3d;
425
426   for ( ii = 0 ; ii < NCONTROL; ii++) {
427     pcons[ii] = wcons;
428     pc3d[ii]  = wc3d;
429     wcons += deltacons;
430     wc3d  += deltac3d;
431   }
432   pcons[NCONTROL] = lcons;
433   pc3d[NCONTROL]  = lc3d;
434
435   Standard_Integer New_NCONTROL = NCONTROL;
436   if(Continuity < GeomAbs_C1) {
437      Standard_Integer NbInt = myHCurve2d->NbIntervals(GeomAbs_C1) + 1;
438      TColStd_Array1OfReal Param_de_decoupeC1 (1, NbInt);
439      myHCurve2d->Intervals(Param_de_decoupeC1, GeomAbs_C1);
440      TColStd_SequenceOfReal new_par;
441      Standard_Integer inter = 1;
442      ii =1;
443      new_par.Append(fcons);
444
445      while(Param_de_decoupeC1(inter) <= fcons + deltamin) inter++;
446      while(Param_de_decoupeC1(NbInt) >= lcons - deltamin) NbInt--;
447
448      while(inter <= NbInt || ii < NCONTROL) {
449        if(Param_de_decoupeC1(inter) < pcons[ii]) {
450          new_par.Append(Param_de_decoupeC1(inter));
451          if((pcons[ii] - Param_de_decoupeC1(inter)) <= deltamin) {
452            ii++;
453            if(ii > NCONTROL) {ii = NCONTROL;}
454          }
455          inter++;
456        }
457        else {
458          if((Param_de_decoupeC1(inter) - pcons[ii]) > deltamin) {
459            new_par.Append(pcons[ii]);
460          }
461          ii++;
462        }
463      }
464      
465      new_par.Append(lcons);
466      New_NCONTROL = new_par.Length() - 1;
467      //simple protection if New_NCONTROL > allocated elements in array
468      if (New_NCONTROL > MAX_ARRAY_SIZE) {
469        mySameParameter = Standard_False;
470        return;
471      }
472      for(ii = 1; ii <= New_NCONTROL; ii++){
473        pcons[ii] = pc3d[ii] = new_par.Value(ii + 1);
474      }
475      pc3d[New_NCONTROL]  = lc3d;
476    }
477
478   
479   Extrema_LocateExtPC Projector;
480   Projector.Initialize(myC3d->Curve(),fc3d,lc3d,Tol);
481   
482   Standard_Integer count = 1;
483   Standard_Real previousp = fc3d, initp=0, curp;//, deltamin = 50*Tolp;
484   Standard_Real bornesup = lc3d - deltamin;
485   Standard_Boolean projok = 0, 
486     use_parameter ;
487   for (ii = 1; ii < New_NCONTROL; ii++){    
488     CurveOnSurface.D0(pcons[ii],Pcons);
489     myC3d->D0(pc3d[ii],Pc3d);
490     dist2 = Pcons.SquareDistance(Pc3d);
491     use_parameter = (dist2 <= Tol2  && (pc3d[ii] > pc3d[count-1] + deltamin)) ;
492     if(use_parameter) {
493       
494       if(dist2 > dmax2) dmax2 = dist2;
495       initp = previousp = pc3d[count] = pc3d[ii];
496       pcons[count] = pcons[ii];
497       count++;
498     }
499     else {
500       if(!projok) initp = pc3d[ii];
501       projok = mySameParameter = Standard_False;
502       Projector.Perform(Pcons, initp);
503       if (Projector.IsDone()) {
504         curp = Projector.Point().Parameter();
505         Standard_Real dist_2 = Projector.SquareDistance();
506         if(dist_2 > besttol2) besttol2 = dist_2;
507         projok = 1;
508       }
509       else {
510         ProjectPointOnCurve(initp,Pcons,Tol,30,myC3d->Curve(),projok,curp);
511       }
512       if(projok){
513         if(curp > previousp + deltamin && curp < bornesup){
514           initp = previousp = pc3d[count] = curp;
515           pcons[count] = pcons[ii];
516           count++;
517         }
518       }
519       else {
520 #ifdef DEB 
521         // JAG
522         cout << "Projection not done" << endl;
523 #endif
524       }
525     }
526   }
527   if(mySameParameter){
528     myTolReached = 1.5*sqrt(dmax2);
529     return;
530   }
531  
532   if(!extrok) { // If not already SameP and tangent to mill, abandon.
533     mySameParameter = Standard_False;
534 #ifdef DEB
535     cout<<"SameParameter problem  : zero tangent to extremities"<<endl;
536 #endif
537     return;
538   }
539
540   pcons[count] = lcons;
541   pc3d[count]  = lc3d;
542
543 #ifdef DEB
544   if (AffichFw) {
545     char Name[17];
546     Name[0]='\0';
547     TColgp_Array1OfPnt2d    DEBP2d  (0,count);
548     TColStd_Array1OfInteger DEBMults(0,count); 
549     DEBMults.Init(1); DEBMults(0) = 2; DEBMults(count) = 2;
550     TColStd_Array1OfReal    DEBKnots(0,count);
551     for (Standard_Integer DEBi = 0; DEBi <= count; DEBi++) {
552       DEBKnots(DEBi) = DEBi;
553       DEBP2d  (DEBi) = gp_Pnt2d(pc3d[DEBi],pcons[DEBi]);
554     }
555     Handle(Geom2d_BSplineCurve) DEBBS = 
556       new Geom2d_BSplineCurve(DEBP2d,DEBKnots,DEBMults,1);
557     sprintf(Name,"DEBC2d_%d",++NbCurve);
558 #ifdef DRAW
559     DrawTrSurf::Set(Name,DEBBS);
560 #endif
561   }
562 #endif
563     
564   while(!interpolok){
565     // The tables and their limits for the interpolation.
566     Standard_Integer num_knots = count + 7;
567     Standard_Integer num_poles = count + 3;
568     TColStd_Array1OfReal    Paramc3d(*pc3d,1,count+1);
569     TColStd_Array1OfReal    Paramcons(*pcons,1,count+1);
570     TColStd_Array1OfInteger ContactOrder(1,num_poles) ;
571     TColStd_Array1OfReal    Poles(1,num_poles) ;
572     TColStd_Array1OfReal    InterpolationParameters(1,num_poles) ;
573     TColStd_Array1OfReal    FlatKnots(1,num_knots) ; 
574     
575     // Fill tables taking attention to end values.
576     ContactOrder.Init(0);
577     ContactOrder(2) = ContactOrder(num_poles - 1) = 1;
578     
579     FlatKnots(1) = FlatKnots(2) = FlatKnots(3) = FlatKnots(4) = fc3d;
580     FlatKnots(num_poles + 1) = FlatKnots(num_poles + 2) = 
581       FlatKnots(num_poles + 3) = FlatKnots(num_poles + 4) = lc3d;
582     
583     Poles(1) = fcons; Poles(num_poles) = lcons;
584     Poles(2) = tangent[0]; Poles(num_poles - 1) = tangent[1];
585     
586     InterpolationParameters(1) = InterpolationParameters(2) = fc3d;
587     InterpolationParameters(num_poles - 1) = InterpolationParameters(num_poles) = lc3d;
588     
589     for (ii = 3; ii <= num_poles - 2; ii++) {
590       Poles(ii) = Paramcons(ii - 1);
591       InterpolationParameters(ii) = FlatKnots(ii+2) = Paramc3d(ii - 1);
592     }
593     Standard_Integer inversion_problem;
594     BSplCLib::Interpolate(3,FlatKnots,InterpolationParameters,ContactOrder,
595                           1,Poles(1),inversion_problem);
596     if(inversion_problem) {
597       Standard_ConstructionError::Raise();
598     }
599
600     //-------------------------------------------
601 #ifdef DEB
602   if (AffichFw) {
603     nbcoups ++;
604     char Name[17];
605     Name[0] = '\0';
606     Standard_Integer nnn = 100;
607     TColgp_Array1OfPnt2d    DEBP2d  (0,nnn);
608     TColStd_Array1OfInteger DEBMults(0,nnn); 
609     DEBMults.Init(1); DEBMults(0) = 2; DEBMults(nnn) = 2;
610     TColStd_Array1OfReal    DEBKnots(0,nnn);
611     Standard_Real du = (lc3d - fc3d) / nnn;
612     Standard_Real u3d = fc3d;
613     Standard_Integer extrap_mode[2] ;
614     extrap_mode[0] = extrap_mode[1] = 3;
615     Standard_Real eval_result[2] ;
616     Standard_Integer DerivativeRequest = 0;
617     Standard_Real *PolesArray =
618       (Standard_Real *) &Poles(Poles.Lower()) ;
619
620     for (Standard_Integer DEBi = 0; DEBi <= nnn; DEBi++) {
621       DEBKnots(DEBi) = DEBi;
622       BSplCLib::Eval(u3d,
623                      Standard_False,
624                      DerivativeRequest,
625                      extrap_mode[0],
626                      3,
627                      FlatKnots,
628                      1,
629                      PolesArray[0],
630                      eval_result[0]) ;
631
632       DEBP2d  (DEBi) = gp_Pnt2d(u3d,eval_result[0]);
633       u3d += du;
634     }
635
636     Handle(Geom2d_BSplineCurve) DEBBS = 
637       new Geom2d_BSplineCurve(DEBP2d,DEBKnots,DEBMults,1);
638     sprintf(Name,"DEBC2d_%d_%d",NbCurve,nbcoups );
639 #ifdef DRAW
640     DrawTrSurf::Set(Name,DEBBS);
641 #endif
642   }
643 #endif
644 //-------------------------------------------    
645
646 //-------------------------------------------
647     // Test if par2d(par3d) is monotonous function or not           ----- IFV, Jan 2000
648     // and try to insert new point to improve BSpline interpolation
649
650     Standard_Integer extrap_mode[2] ;
651     extrap_mode[0] = extrap_mode[1] = 3;
652     Standard_Real eval_result[2] ;
653     Standard_Integer DerivativeRequest = 0;
654     Standard_Real *PolesArray =
655       (Standard_Real *) &Poles(Poles.Lower()) ;
656
657     Standard_Integer newcount = 0;
658     for (ii = 0; ii < count; ii++) {
659       
660       newpcons[newcount] = pcons[ii];
661       newpc3d[newcount] = pc3d[ii];
662       newcount++;
663
664       if(count - ii + newcount == MAX_ARRAY_SIZE) continue;
665
666       BSplCLib::Eval(.5*(pc3d[ii]+pc3d[ii+1]), Standard_False, DerivativeRequest,
667                      extrap_mode[0], 3, FlatKnots, 1, PolesArray[0], eval_result[0]);
668                      
669       if(eval_result[0] < pcons[ii] || eval_result[0] > pcons[ii+1]) {
670         Standard_Real ucons = 0.5*(pcons[ii]+pcons[ii+1]);
671         Standard_Real uc3d  = 0.5*(pc3d[ii]+pc3d[ii+1]);
672         
673         CurveOnSurface.D0(ucons,Pcons);
674         Projector.Perform(Pcons, uc3d);
675         if (Projector.IsDone()) {
676           curp = Projector.Point().Parameter();
677           Standard_Real dist_2 = Projector.SquareDistance();
678           if(dist_2 > besttol2) besttol2 = dist_2;
679           projok = 1;
680         }
681         else {
682           ProjectPointOnCurve(uc3d,Pcons,Tol,30,myC3d->Curve(),projok,curp);
683         }
684         if(projok){
685           if(curp > pc3d[ii] + deltamin && curp < pc3d[ii+1] - deltamin){
686             newpc3d[newcount] = curp;
687             newpcons[newcount] = ucons;
688             newcount ++;
689           }
690         }
691         else {
692 #ifdef DEB 
693           // JAG
694           cout << "Projection not done" << endl;
695 #endif
696         }
697       }
698      
699     }
700
701     newpc3d[newcount] = pc3d[count];
702     newpcons[newcount] = pcons[count];
703     Standard_Real * temp;
704     temp = pc3d;
705     pc3d = newpc3d;
706     newpc3d = temp;
707     temp = pcons;
708     pcons = newpcons;
709     newpcons = temp;
710
711     if((count != newcount) && newcount < MAX_ARRAY_SIZE) { count = newcount; continue;}
712
713     count = newcount;
714
715     Standard_Real algtol = sqrt(besttol2);
716
717     interpolok = Check (FlatKnots, Poles, count+1, Paramc3d, Paramcons, 
718                         myC3d, CurveOnSurface, algtol, tolsov);
719
720     if (Precision::IsInfinite(algtol)) {
721       mySameParameter = Standard_False;
722 #ifdef DEB
723       cout<<"SameParameter problem  : function of interpolation of parametration at mills !!"<<endl;
724 #endif
725       return;
726     }
727
728     tolsov = algtol;
729
730     interpolok = (interpolok || count >= MAX_ARRAY_SIZE);
731
732     if(interpolok) {
733         Standard_Real besttol = sqrt(besttol2);
734 #ifdef DEB
735       if (Voir) {
736         if(algtol > besttol){
737           cout<<"SameParameter : Tol can't be reached before approx"<<endl;
738         }
739       }
740 #endif
741       Handle(TColStd_HArray1OfReal) tol1d,tol2d,tol3d;
742       tol1d = new TColStd_HArray1OfReal(1,2) ;
743       tol1d->SetValue(1, mySurf->UResolution(besttol));
744       tol1d->SetValue(2, mySurf->VResolution(besttol));
745
746       Approx_SameParameter_Evaluator ev (FlatKnots, Poles, myHCurve2d); 
747       AdvApprox_ApproxAFunction  anApproximator(2,0,0,tol1d,tol2d,tol3d,fc3d,lc3d,
748                                                 Continuity,11,40,ev);
749
750       if (anApproximator.IsDone() || anApproximator.HasResult()) {
751         GeomLib_MakeCurvefromApprox  aCurveBuilder(anApproximator) ;
752         myCurve2d = aCurveBuilder.Curve2dFromTwo1d(1,2) ;
753         myHCurve2d = new Geom2dAdaptor_HCurve(myCurve2d);
754         CurveOnSurface.Load(myHCurve2d);
755         myTolReached = ComputeTolReached(myC3d,CurveOnSurface,NCONTROL);
756         myDone = Standard_True;
757       } 
758     }
759     else {
760 #ifdef DEB
761       if (Voir)
762         cout<<"SameParameter : Not enough points, enrich"<<endl;
763 #endif
764
765       Standard_Integer newcount = 0;
766       for(Standard_Integer n = 0; n < count; n++){
767         newpc3d[newcount] = pc3d[n];
768         newpcons[newcount] = pcons[n];
769         newcount ++;
770
771         if(count - n + newcount == MAX_ARRAY_SIZE) continue;
772
773         Standard_Real ucons = 0.5*(pcons[n]+pcons[n+1]);
774         Standard_Real uc3d  = 0.5*(pc3d[n]+pc3d[n+1]);
775         
776         CurveOnSurface.D0(ucons,Pcons);
777         Projector.Perform(Pcons, uc3d);
778         if (Projector.IsDone()) {
779           curp = Projector.Point().Parameter();
780           Standard_Real dist_2 = Projector.SquareDistance();
781           if(dist_2 > besttol2) besttol2 = dist_2;
782           projok = 1;
783         }
784         else {
785           ProjectPointOnCurve(uc3d,Pcons,Tol,30,myC3d->Curve(),projok,curp);
786         }
787         if(projok){
788           if(curp > pc3d[n] + deltamin && curp < pc3d[n+1] - deltamin){
789             newpc3d[newcount] = curp;
790             newpcons[newcount] = ucons;
791             newcount ++;
792           }
793         }
794         else {
795 #ifdef DEB 
796           // JAG
797           cout << "Projection not done" << endl;
798 #endif
799         }
800       }
801       newpc3d[newcount] = pc3d[count];
802       newpcons[newcount] = pcons[count];
803       Standard_Real * tempx;
804       tempx = pc3d;
805       pc3d = newpc3d;
806       newpc3d = tempx;
807       tempx = pcons;
808       pcons = newpcons;
809       newpcons = tempx;
810       count = newcount;
811     }
812   }
813 }