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