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