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