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