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