7258312053da118a52dcb1ceb2ea76e69b0e1850
[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
18 #include <Adaptor2d_HCurve2d.hxx>
19 #include <Adaptor3d_CurveOnSurface.hxx>
20 #include <Adaptor3d_HCurve.hxx>
21 #include <Adaptor3d_HSurface.hxx>
22 #include <AdvApprox_ApproxAFunction.hxx>
23 #include <Approx_SameParameter.hxx>
24 #include <BSplCLib.hxx>
25 #include <Extrema_ExtPC.hxx>
26 #include <Extrema_LocateExtPC.hxx>
27 #include <Geom2d_BSplineCurve.hxx>
28 #include <Geom2d_Curve.hxx>
29 #include <Geom2dAdaptor_Curve.hxx>
30 #include <Geom2dAdaptor_HCurve.hxx>
31 #include <Geom_Curve.hxx>
32 #include <Geom_Surface.hxx>
33 #include <GeomAdaptor_HCurve.hxx>
34 #include <GeomAdaptor_HSurface.hxx>
35 #include <GeomLib_MakeCurvefromApprox.hxx>
36 #include <Precision.hxx>
37 #include <Standard_ConstructionError.hxx>
38 #include <Standard_OutOfRange.hxx>
39 #include <TColStd_Array1OfReal.hxx>
40 #include <Geom2dAdaptor.hxx>
41
42 //=======================================================================
43 //class    : Approx_SameParameter_Evaluator
44 //purpose  : Used in same parameterization curve approximation.
45 //=======================================================================
46 class Approx_SameParameter_Evaluator : public AdvApprox_EvaluatorFunction
47 {
48 public:
49   Approx_SameParameter_Evaluator (const TColStd_Array1OfReal& theFlatKnots,
50                                   const TColStd_Array1OfReal& thePoles,
51                                   const Handle(Adaptor2d_HCurve2d)& theHCurve2d)
52     : FlatKnots(theFlatKnots),
53       Poles(thePoles),
54       HCurve2d(theHCurve2d) {}
55
56   virtual void Evaluate (Standard_Integer *Dimension,
57                          Standard_Real     StartEnd[2],
58                          Standard_Real    *Parameter,
59                          Standard_Integer *DerivativeRequest,
60                          Standard_Real    *Result, // [Dimension]
61                          Standard_Integer *ErrorCode);
62
63 private:
64   const TColStd_Array1OfReal& FlatKnots;
65   const TColStd_Array1OfReal& Poles;
66   Handle(Adaptor2d_HCurve2d) HCurve2d;
67 };
68
69 //=======================================================================
70 //function : Evaluate
71 //purpose  : 
72 //=======================================================================
73 void Approx_SameParameter_Evaluator::Evaluate (Standard_Integer *,/*Dimension*/
74                                                Standard_Real    /*StartEnd*/[2],
75                                                Standard_Real    *Parameter,
76                                                Standard_Integer *DerivativeRequest,
77                                                Standard_Real    *Result,
78                                                Standard_Integer *ReturnCode) 
79 {
80   const Standard_Integer aDegree = 3;
81   Standard_Integer extrap_mode[2] = {aDegree, aDegree};
82   Standard_Real eval_result[2];
83   Standard_Real *PolesArray = (Standard_Real *) &Poles(Poles.Lower()) ;
84
85   // Evaluate the 1D B-Spline that represents the change in parameterization.
86   BSplCLib::Eval(*Parameter,
87                  Standard_False,
88                  *DerivativeRequest,
89                  extrap_mode[0],
90                  aDegree,
91                  FlatKnots,
92                  1,
93                  PolesArray[0],
94                  eval_result[0]);
95
96   gp_Pnt2d aPoint;
97   gp_Vec2d aVector;
98   if (*DerivativeRequest == 0)
99   {
100     HCurve2d->D0(eval_result[0], aPoint);
101     aPoint.Coord(Result[0],Result[1]);
102   }
103   else if (*DerivativeRequest == 1)
104   {
105     HCurve2d->D1(eval_result[0], aPoint, aVector);
106     aVector.Multiply(eval_result[1]);
107     aVector.Coord(Result[0],Result[1]);
108   }
109
110   ReturnCode[0] = 0;
111 }
112
113 //=======================================================================
114 //function : ProjectPointOnCurve
115 //purpose  : 
116 //=======================================================================
117 static void ProjectPointOnCurve(const Standard_Real      InitValue,
118                                 const gp_Pnt             APoint,
119                                 const Standard_Real      Tolerance,
120                                 const Standard_Integer   NumIteration,
121                                 const Adaptor3d_Curve&   Curve,
122                                 Standard_Boolean&        Status,
123                                 Standard_Real&           Result)
124 {
125   Standard_Integer num_iter = 0, not_done = 1;
126
127   gp_Pnt a_point;
128   gp_Vec vector, d1, d2;
129   Standard_Real func, func_derivative,
130                 param = InitValue;
131   Status = Standard_False;
132   do
133   {
134     num_iter++;
135     Curve.D2(param, a_point, d1, d2);
136     vector = gp_Vec(a_point,APoint);
137
138     func = vector.Dot(d1);
139     if ( Abs(func) < Tolerance * d1.Magnitude())
140     {
141       not_done = 0;
142       Status = Standard_True;
143     }
144     else
145     {
146       func_derivative = vector.Dot(d2) - d1.Dot(d1);
147
148       // Avoid division by zero.
149       const Standard_Real Toler = 1.0e-12;
150       if( Abs(func_derivative) > Toler )
151         param -= func / func_derivative;
152
153       param = Max(param,Curve.FirstParameter());
154       param = Min(param,Curve.LastParameter());
155     }
156   } while (not_done && num_iter <= NumIteration);
157
158   Result = param;
159 }
160
161 //=======================================================================
162 //function : ComputeTolReached
163 //purpose  :
164 //=======================================================================
165 static Standard_Real ComputeTolReached(const Handle(Adaptor3d_HCurve)& c3d,
166                                        const Adaptor3d_CurveOnSurface& cons,
167                                        const Standard_Integer        nbp)
168 {
169   Standard_Real d2 = 0.0; // Square max discrete deviation.
170   const Standard_Real first = c3d->FirstParameter();
171   const Standard_Real last  = c3d->LastParameter();
172   for(Standard_Integer i = 0; i <= nbp; i++)
173   {
174     Standard_Real t = IntToReal(i) / IntToReal(nbp);
175     Standard_Real u = first * (1.0 - t) + last * t;
176     gp_Pnt Pc3d = c3d->Value(u);
177     gp_Pnt Pcons = cons.Value(u);
178     if (Precision::IsInfinite(Pcons.X()) ||
179         Precision::IsInfinite(Pcons.Y()) ||
180         Precision::IsInfinite(Pcons.Z()))
181     {
182         d2=Precision::Infinite();
183         break;
184     }
185     d2 = Max(d2, Pc3d.SquareDistance(Pcons));
186   }
187
188   const Standard_Real aMult = 1. + 0.05; // 
189   Standard_Real aDeviation = aMult * sqrt(d2);
190   aDeviation = Max(aDeviation, Precision::Confusion()); // Tolerance in modeling space.
191   return aDeviation;
192 }
193
194 //=======================================================================
195 //function : Check
196 //purpose  : Check current interpolation for validity.
197 //=======================================================================
198 static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots,
199                               const TColStd_Array1OfReal& Poles,
200                               const Standard_Integer nbp,
201                               const Standard_Real *pc3d,
202                               const Handle(Adaptor3d_HCurve)& c3d,
203                               const Adaptor3d_CurveOnSurface& cons,
204                               Standard_Real& tol,
205                               const Standard_Real oldtol)
206 {
207   const Standard_Integer aDegree = 3;
208   Standard_Integer extrap_mode[2] = {aDegree, aDegree};
209
210   // Correction of the interval of valid values. This condition has no sensible
211   // grounds. But it is better then the old one (which is commented out) because
212   // it fixes the bug OCC5898. To develop more or less sensible criterion it is
213   // necessary to deeply investigate this problem which is not possible in frames
214   // of debugging.
215   Standard_Real aParamFirst = 3.0 * pc3d[0]   - 2.0 * pc3d[nbp - 1];
216   Standard_Real aParamLast = 3.0 * pc3d[nbp - 1] - 2.0 * pc3d[0];
217
218   Standard_Real FirstPar = cons.FirstParameter();
219   Standard_Real LastPar  = cons.LastParameter();
220   if (aParamFirst < FirstPar)
221     aParamFirst = FirstPar;
222   if (aParamLast > LastPar)
223     aParamLast = LastPar;
224
225   Standard_Real d2 = 0.0; // Maximum square deviation on the samples.
226   const Standard_Real d = tol;
227   const Standard_Integer nn = 2 * nbp;
228   const Standard_Real unsurnn = 1.0 / nn;
229   Standard_Real tprev = aParamFirst;
230   for(Standard_Integer i = 0; i <= nn; i++)
231   {
232     // Compute corresponding parameter on 2d curve.
233     // It should be inside of 3d curve parameter space.
234     Standard_Real t = unsurnn*i;
235     Standard_Real tc3d = pc3d[0]*(1.0 - t) + pc3d[nbp - 1] * t; // weight function.
236     gp_Pnt Pc3d = c3d->Value(tc3d);
237     Standard_Real tcons;
238     BSplCLib::Eval(tc3d, Standard_False, 0, extrap_mode[0],
239                    aDegree, FlatKnots, 1, (Standard_Real&)Poles(1), tcons);
240
241     if (tcons < tprev ||
242         tcons > aParamLast)
243     {
244       tol = Precision::Infinite();
245       return Standard_False;
246     }
247     tprev = tcons;
248     gp_Pnt Pcons = cons.Value(tcons);
249     Standard_Real temp = Pc3d.SquareDistance(Pcons);
250     if(temp > d2) d2 = temp;
251   }
252   tol = sqrt(d2);
253
254   // Check poles parameters to be ordered.
255   for(Standard_Integer i = Poles.Lower() + 1; i <= Poles.Upper(); ++i)
256   {
257     const Standard_Real aPreviousParam = Poles(i - 1);
258     const Standard_Real aCurrentParam  = Poles(i);
259
260     if (aPreviousParam > aCurrentParam)
261       return Standard_False;
262   }
263
264   return (tol <= d || tol > 0.8 * oldtol);
265 }
266
267 //=======================================================================
268 //function : Approx_SameParameter
269 //purpose  : 
270 //=======================================================================
271 Approx_SameParameter::Approx_SameParameter(const Handle(Geom_Curve)&   C3D,
272                                            const Handle(Geom2d_Curve)& C2D,
273                                            const Handle(Geom_Surface)& S,
274                                            const Standard_Real         Tol)
275 : myDeltaMin(Precision::PConfusion()),
276   mySameParameter(Standard_True),
277   myDone(Standard_False)
278 {
279   myHCurve2d = new Geom2dAdaptor_HCurve(C2D);
280   myC3d      = new GeomAdaptor_HCurve(C3D);
281   mySurf     = new GeomAdaptor_HSurface(S);
282   Build(Tol);
283 }
284
285 //=======================================================================
286 //function : Approx_SameParameter
287 //purpose  : 
288 //=======================================================================
289 Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)&   C3D,
290                                            const Handle(Geom2d_Curve)&       C2D,
291                                            const Handle(Adaptor3d_HSurface)& S,
292                                            const Standard_Real               Tol)
293 : myDeltaMin(Precision::PConfusion()),
294   mySameParameter(Standard_True),
295   myDone(Standard_False)
296 {
297   myC3d = C3D;
298   mySurf = S;
299   myHCurve2d = new Geom2dAdaptor_HCurve(C2D);
300   Build(Tol);
301 }
302
303 //=======================================================================
304 //function : Approx_SameParameter
305 //purpose  : 
306 //=======================================================================
307 Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)&   C3D,
308                                            const Handle(Adaptor2d_HCurve2d)& C2D,
309                                            const Handle(Adaptor3d_HSurface)& S,
310                                            const Standard_Real               Tol)
311 : myDeltaMin(Precision::PConfusion()),
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   // Algorithm:
328   // 1) Build initial distribution. Increase number of samples in case of C0 continuity.
329   // 2.1) Check same parameter state on samples.
330   // 2.2) Compute parameters in 2d space if not same parameter.
331   // 3) Loop over poles number and try to interpolate 2d curve.
332   // 4) If loop is failed build 2d curve forcibly or use original pcurve.
333   Standard_Real qpcons[myMaxArraySize], qnewpcons[myMaxArraySize],
334                  qpc3d[myMaxArraySize], qnewpc3d[myMaxArraySize];
335
336   // Create and fill data structure.
337   Approx_SameParameter_Data aData;
338   aData.myCOnS = Adaptor3d_CurveOnSurface(myHCurve2d,mySurf);
339   aData.myC2dPF = aData.myCOnS.FirstParameter();
340   aData.myC2dPL = aData.myCOnS.LastParameter();
341   aData.myC3dPF = myC3d->FirstParameter();
342   aData.myC3dPL = myC3d->LastParameter();
343   aData.myNbPnt = 0; // No points initially.
344   aData.myPC2d = qpcons;
345   aData.myPC3d = qpc3d;
346   aData.myNewPC2d = qnewpcons;
347   aData.myNewPC3d = qnewpc3d;
348   aData.myTol = Tolerance;
349
350   // Build initial distribution.
351   if (!BuildInitialDistribution(aData))
352   {
353     mySameParameter = Standard_False;
354     myDone = Standard_False;
355     return;
356   }
357
358   // Check same parameter state on distribution.
359   Standard_Real aMaxSqDeviation = 0.0;
360   const Standard_Real aPercentOfBadProj = 0.3;
361   Standard_Integer aNbPnt = aData.myNbPnt - RealToInt(aPercentOfBadProj * aData.myNbPnt);
362   mySameParameter = CheckSameParameter(aData, aMaxSqDeviation);
363   if(mySameParameter)
364   {
365     myTolReached = ComputeTolReached(myC3d, aData.myCOnS, 2 * myNbSamples);
366     myDone = Standard_True;
367     return;
368   }
369   else
370   {
371     // Control number of sample points after checking sameparameter
372     // If number of points is less then initial one, it means that there are 
373     // problems with projection
374     if(aData.myNbPnt < aNbPnt )
375     {
376       myTolReached = ComputeTolReached(myC3d,aData.myCOnS, 2 * myNbSamples);
377       myCurve2d = Geom2dAdaptor::MakeCurve(myHCurve2d->Curve2d());
378       myDone = Standard_False;
379       return;
380     }
381   }
382
383   // Control tangents at the extremities to know if the
384   // reparametring is possible and calculate the tangents
385   // at the extremities of the function of change of variable.
386   Standard_Real tangent[2] = { 0.0, 0.0 };
387   if(!ComputeTangents(aData.myCOnS, tangent[0], tangent[1]))
388   {
389     // Cannot compute tangents.
390     mySameParameter = Standard_False;
391     myDone = Standard_False;
392     myTolReached = ComputeTolReached(myC3d, aData.myCOnS, 2 * myNbSamples);
393     return;
394   }
395
396   // There is at least one point where same parameter is broken.
397   // Try to build B-spline approximation curve using interpolation with degree 3.
398   // The loop is organized over number of poles.
399   GeomAbs_Shape aContinuity = myHCurve2d->Continuity();
400   if(aContinuity > GeomAbs_C1) aContinuity = GeomAbs_C1;
401
402   Standard_Real besttol2 = aData.myTol * aData.myTol,
403                 tolsov = Precision::Infinite();
404   Standard_Boolean interpolok = Standard_False,
405                    hasCountChanged = Standard_False;
406   do
407   {
408     // Interpolation data.
409     Standard_Integer num_knots = aData.myNbPnt + 7;
410     Standard_Integer num_poles = aData.myNbPnt + 3;
411     TColStd_Array1OfReal    Poles(1, num_poles);
412     TColStd_Array1OfReal    FlatKnots(1 ,num_knots);
413
414     if (!Interpolate(aData, tangent[0], tangent[1], Poles, FlatKnots))
415     {
416       // Interpolation fails.
417       mySameParameter = Standard_False;
418       myDone = Standard_False;
419       myTolReached = ComputeTolReached(myC3d, aData.myCOnS, 2 * myNbSamples);
420       return;
421     }
422
423     Standard_Real algtol = sqrt(besttol2);
424     interpolok = Check (FlatKnots, Poles, aData.myNbPnt+1, aData.myPC3d,
425                         myC3d, aData.myCOnS, algtol, tolsov);
426     tolsov = algtol;
427
428     // Try to build 2d curve and check it for validity.
429     if(interpolok)
430     {
431       Standard_Real besttol = sqrt(besttol2);
432
433       Handle(TColStd_HArray1OfReal) tol1d,tol2d,tol3d;
434       tol1d = new TColStd_HArray1OfReal(1,2);
435       tol1d->SetValue(1, mySurf->UResolution(besttol));
436       tol1d->SetValue(2, mySurf->VResolution(besttol));
437
438       Approx_SameParameter_Evaluator ev (FlatKnots, Poles, myHCurve2d);
439       Standard_Integer aMaxDeg = 11, aMaxSeg = 1000;
440       AdvApprox_ApproxAFunction  anApproximator(2,0,0,tol1d,tol2d,tol3d,aData.myC3dPF,aData.myC3dPL,
441                                                 aContinuity,aMaxDeg,aMaxSeg,ev);
442
443       if (anApproximator.IsDone() || anApproximator.HasResult())
444       {
445         Adaptor3d_CurveOnSurface ACS = aData.myCOnS;
446         GeomLib_MakeCurvefromApprox  aCurveBuilder(anApproximator);
447         Handle(Geom2d_BSplineCurve) aC2d = aCurveBuilder.Curve2dFromTwo1d(1,2);
448         Handle(Adaptor2d_HCurve2d) aHCurve2d = new Geom2dAdaptor_HCurve(aC2d);
449         aData.myCOnS.Load(aHCurve2d);
450         myTolReached = ComputeTolReached(myC3d,aData.myCOnS, 2 * myNbSamples);
451
452         const Standard_Real aMult = 250.0; // To be tolerant with discrete tolerance.
453         if (myTolReached < aMult * besttol ) 
454         {
455           myCurve2d = aC2d;
456           myHCurve2d = aHCurve2d;
457           myDone = Standard_True;
458           break;
459         }
460         else if(aData.myNbPnt < myMaxArraySize - 1)
461         {
462           interpolok = Standard_False;
463           aData.myCOnS = ACS;
464         }
465         else
466         {
467           break;
468         }
469       }
470     }
471
472     if (!interpolok)
473       hasCountChanged = IncreaseNbPoles(Poles, FlatKnots, aData, besttol2);
474
475   } while(!interpolok && hasCountChanged);
476
477   if (!myDone)
478   {
479     // Loop is finished unsuccessfully. Fix tolerance by maximal deviation,
480     // using data from the last loop iteration or initial data. Use data set with minimal deflection.
481
482     // Original 2d curve.
483     aData.myCOnS.Load(myHCurve2d);
484     myTolReached = ComputeTolReached(myC3d,aData.myCOnS, 2 * myNbSamples);
485     myCurve2d = Geom2dAdaptor::MakeCurve(myHCurve2d->Curve2d());
486
487     // Approximation curve.
488     Standard_Integer num_knots = aData.myNbPnt + 7;
489     Standard_Integer num_poles = aData.myNbPnt + 3;
490     TColStd_Array1OfReal    Poles(1, num_poles);
491     TColStd_Array1OfReal    FlatKnots(1 ,num_knots);
492
493     Interpolate(aData, tangent[0], tangent[1],
494                 Poles, FlatKnots);
495
496     Standard_Real besttol = sqrt(besttol2);
497     Handle(TColStd_HArray1OfReal) tol1d,tol2d,tol3d;
498     tol1d = new TColStd_HArray1OfReal(1,2) ;
499     tol1d->SetValue(1, mySurf->UResolution(besttol));
500     tol1d->SetValue(2, mySurf->VResolution(besttol));
501
502     Approx_SameParameter_Evaluator ev(FlatKnots, Poles, myHCurve2d);
503     AdvApprox_ApproxAFunction  anApproximator(2,0,0,tol1d,tol2d,tol3d,aData.myC3dPF,aData.myC3dPL,
504                                               aContinuity,11,40,ev);
505
506     if (!anApproximator.IsDone() &&
507         !anApproximator.HasResult() )
508     {
509       myDone = Standard_False;
510       return;
511     }
512
513     GeomLib_MakeCurvefromApprox  aCurveBuilder(anApproximator);
514     Handle(Geom2d_BSplineCurve) aC2d = aCurveBuilder.Curve2dFromTwo1d(1,2);
515     Handle(Adaptor2d_HCurve2d) aHCurve2d = new Geom2dAdaptor_HCurve(aC2d);
516     aData.myCOnS.Load(aHCurve2d);
517
518     Standard_Real anApproxTol = ComputeTolReached(myC3d,aData.myCOnS,2 * myNbSamples);
519     if (anApproxTol < myTolReached)
520     {
521       myTolReached = anApproxTol;
522       myCurve2d = aC2d;
523       myHCurve2d  = aHCurve2d;
524     }
525     myDone = Standard_True;
526   }
527 }
528
529 //=======================================================================
530 //function : BuildInitialDistribution
531 //purpose  : Sub-method in Build.
532 //=======================================================================
533 Standard_Boolean Approx_SameParameter::BuildInitialDistribution(Approx_SameParameter_Data &theData) const
534 {
535   // Take a multiple of the sample pof CheckShape,
536   // at least the control points will be correct.
537   // Take parameters with constant step on the curve on surface
538   // and on curve 3d.
539   const Standard_Real deltacons = (theData.myC2dPL - theData.myC2dPF) / myNbSamples;
540   const Standard_Real deltac3d  = (theData.myC3dPL - theData.myC3dPF) / myNbSamples;
541   Standard_Real wcons = theData.myC2dPF;
542   Standard_Real wc3d  = theData.myC3dPF;
543   for (Standard_Integer ii = 0 ; ii < myNbSamples; ii++)
544   {
545     theData.myPC2d[ii] = wcons;
546     theData.myPC3d[ii]  = wc3d;
547     wcons += deltacons;
548     wc3d  += deltac3d;
549   }
550   theData.myNbPnt = myNbSamples;
551   theData.myPC2d[theData.myNbPnt] = theData.myC2dPL;
552   theData.myPC3d[theData.myNbPnt]  = theData.myC3dPL;
553
554   // Change number of points in case of C0 continuity.
555   GeomAbs_Shape Continuity = myHCurve2d->Continuity();
556   if(Continuity < GeomAbs_C1)
557   {
558     if (!IncreaseInitialNbSamples(theData))
559     {
560       // Number of samples is too big.
561       return Standard_False;
562     }
563   }
564
565   return Standard_True;
566 }
567
568 //=======================================================================
569 //function : IncreaseInitialNbSamples
570 //purpose  : Get number of C1 intervals and build new distribution on them.
571 //           Sub-method in BuildInitialDistribution.
572 //=======================================================================
573 Standard_Boolean Approx_SameParameter::IncreaseInitialNbSamples(Approx_SameParameter_Data &theData) const
574 {
575   Standard_Integer NbInt = myHCurve2d->NbIntervals(GeomAbs_C1) + 1;
576   TColStd_Array1OfReal aC1Intervals (1, NbInt);
577   myHCurve2d->Intervals(aC1Intervals, GeomAbs_C1);
578
579   Standard_Integer inter = 1;
580   while(inter <= NbInt && aC1Intervals(inter) <= theData.myC3dPF + myDeltaMin) inter++;
581   while(NbInt > 0 && aC1Intervals(NbInt) >= theData.myC3dPL - myDeltaMin) NbInt--;
582
583   // Compute new parameters.
584   TColStd_SequenceOfReal aNewPar;
585   aNewPar.Append(theData.myC3dPF);
586   Standard_Integer ii = 1;
587   while(inter <= NbInt || (ii < myNbSamples && inter <= aC1Intervals.Length()) )
588   {
589     if(aC1Intervals(inter) < theData.myPC2d[ii])
590     {
591       aNewPar.Append(aC1Intervals(inter));
592       if((theData.myPC2d[ii] - aC1Intervals(inter)) <= myDeltaMin)
593       {
594         ii++;
595         if(ii > myNbSamples)
596         {
597           ii = myNbSamples;
598         }
599       }
600       inter++;
601     }
602     else
603     {
604       if((aC1Intervals(inter) - theData.myPC2d[ii]) > myDeltaMin)
605       {
606         aNewPar.Append(theData.myPC2d[ii]);
607       }
608       ii++;
609     }
610   }
611   // Simple protection if theNewNbPoints > allocated elements in array but one
612   // myMaxArraySize - 1 index may be filled after projection.
613   theData.myNbPnt = aNewPar.Length();
614   if (theData.myNbPnt > myMaxArraySize - 1)
615   {
616     return Standard_False;
617   }
618
619   for(ii = 1; ii < theData.myNbPnt; ii++)
620   {
621     // Copy only internal points.
622     theData.myPC2d[ii] = theData.myPC3d[ii] = aNewPar.Value(ii + 1);
623   }
624   theData.myPC3d[theData.myNbPnt]  = theData.myC3dPL;
625   theData.myPC2d[theData.myNbPnt] = theData.myC2dPL;
626
627   return Standard_True;
628 }
629
630 //=======================================================================
631 //function : CheckSameParameter
632 //purpose  : Sub-method in Build.
633 //=======================================================================
634 Standard_Boolean Approx_SameParameter::CheckSameParameter(Approx_SameParameter_Data &theData,
635                                                           Standard_Real &theSqDist) const
636 {
637   const Standard_Real Tol2 = theData.myTol * theData.myTol;
638   Standard_Boolean isSameParam = Standard_True;
639
640   // Compute initial distance on boundary points.
641   gp_Pnt Pcons, Pc3d;
642   theData.myCOnS.D0(theData.myC2dPF, Pcons);
643   myC3d->D0(theData.myC3dPF, Pc3d);
644   Standard_Real dist2 = Pcons.SquareDistance(Pc3d);
645   Standard_Real dmax2 = dist2;
646
647   theData.myCOnS.D0(theData.myC2dPL, Pcons);
648   myC3d->D0(theData.myC3dPL, Pc3d);
649   dist2 = Pcons.SquareDistance(Pc3d);
650   dmax2 = Max(dmax2, dist2);
651
652   Extrema_LocateExtPC Projector;
653   Projector.Initialize(myC3d->Curve(), theData.myC3dPF, theData.myC3dPL, theData.myTol);
654
655   Standard_Integer count = 1;
656   Standard_Real previousp = theData.myC3dPF, initp=0, curp;
657   Standard_Real bornesup = theData.myC3dPL - myDeltaMin;
658   Standard_Boolean isProjOk = Standard_False;
659   for (Standard_Integer ii = 1; ii < theData.myNbPnt; ii++)
660   {
661     theData.myCOnS.D0(theData.myPC2d[ii],Pcons);
662     myC3d->D0(theData.myPC3d[ii],Pc3d);
663     dist2 = Pcons.SquareDistance(Pc3d);
664
665     // Same parameter point.
666     Standard_Boolean isUseParam = (dist2 <= Tol2  && // Good distance.
667       (theData.myPC3d[ii] > theData.myPC3d[count-1] + myDeltaMin)); // Point is separated from previous.
668     if(isUseParam)
669     {
670       if(dmax2 < dist2)
671         dmax2 = dist2;
672       initp = previousp = theData.myPC3d[count] = theData.myPC3d[ii];
673       theData.myPC2d[count] = theData.myPC2d[ii];
674       count++;
675       continue;
676     }
677
678     // Local search: local extrema and iterative projection algorithm.
679     if(!isProjOk)
680       initp = theData.myPC3d[ii];
681     isProjOk = isSameParam = Standard_False;
682     Projector.Perform(Pcons, initp);
683     if (Projector.IsDone())
684     {
685       // Local extrema is found.
686       curp = Projector.Point().Parameter();
687       isProjOk = Standard_True;
688     }
689     else
690     {
691       ProjectPointOnCurve(initp,Pcons,theData.myTol,30,myC3d->Curve(),isProjOk,curp);
692     }
693     isProjOk = isProjOk && // Good projection.
694                curp > previousp + myDeltaMin && // Point is separated from previous.
695                curp < bornesup; // Inside of parameter space.
696     if(isProjOk)
697     {
698       initp = previousp = theData.myPC3d[count] = curp;
699       theData.myPC2d[count] = theData.myPC2d[ii];
700       count++;
701       continue;
702     }
703
704     // Whole parameter space search using general extrema.
705     Extrema_ExtPC PR(Pcons,myC3d->Curve(),theData.myC3dPF, theData.myC3dPL,theData.myTol);
706     if (!PR.IsDone() || PR.NbExt() == 0) // Lazy evaluation is used.
707       continue;
708
709     const Standard_Integer aNbExt = PR.NbExt();
710     Standard_Integer anIndMin = 0;
711     Standard_Real aCurDistMin = RealLast();
712     for(Standard_Integer i = 1; i <= aNbExt; i++)
713     {
714       const gp_Pnt &aP = PR.Point(i).Value();
715       Standard_Real aDist2 = aP.SquareDistance(Pcons);
716       if(aDist2 < aCurDistMin)
717       {
718         aCurDistMin = aDist2;
719         anIndMin = i;
720       }
721     }
722     if(anIndMin)
723     {
724       curp = PR.Point(anIndMin).Parameter();
725       if( curp > previousp + myDeltaMin && curp < bornesup)
726       {
727         initp = previousp = theData.myPC3d[count] = curp;
728         theData.myPC2d[count] = theData.myPC2d[ii];
729         count++;
730         isProjOk = Standard_True;
731       }
732     }
733   }
734   theData.myNbPnt = count;
735   theData.myPC2d[theData.myNbPnt] = theData.myC2dPL;
736   theData.myPC3d[theData.myNbPnt] = theData.myC3dPL;
737
738   theSqDist = dmax2;
739   return isSameParam;
740 }
741
742 //=======================================================================
743 //function : ComputeTangents
744 //purpose  : Sub-method in Build.
745 //=======================================================================
746 Standard_Boolean Approx_SameParameter::ComputeTangents(const Adaptor3d_CurveOnSurface & theCOnS,
747                                                        Standard_Real &theFirstTangent,
748                                                        Standard_Real &theLastTangent) const
749 {
750   const Standard_Real aSmallMagnitude = 1.0e-12;
751   // Check tangency on curve border.
752   gp_Pnt aPnt, aPntCOnS;
753   gp_Vec aVec, aVecConS;
754
755   // First point.
756   const Standard_Real aParamFirst = myC3d->FirstParameter();
757   theCOnS.D1(aParamFirst, aPntCOnS, aVecConS);
758   myC3d->D1(aParamFirst, aPnt, aVec);
759   Standard_Real aMagnitude = aVecConS.Magnitude();
760   if (aMagnitude > aSmallMagnitude)
761     theFirstTangent = aVec.Magnitude() / aMagnitude;
762   else
763     return Standard_False;
764
765   // Last point.
766   const Standard_Real aParamLast = myC3d->LastParameter();
767   theCOnS.D1(aParamLast,aPntCOnS,aVecConS);
768   myC3d->D1(aParamLast, aPnt, aVec);
769
770   aMagnitude = aVecConS.Magnitude();
771   if (aMagnitude > aSmallMagnitude)
772     theLastTangent = aVec.Magnitude() / aMagnitude;
773   else
774     return Standard_False;
775
776   return Standard_True;
777 }
778
779 //=======================================================================
780 //function : Interpolate
781 //purpose  : Sub-method in Build.
782 //=======================================================================
783 Standard_Boolean Approx_SameParameter::Interpolate(const Approx_SameParameter_Data & theData,
784                                        const Standard_Real aTangFirst,
785                                        const Standard_Real aTangLast,
786                                        TColStd_Array1OfReal & thePoles,
787                                        TColStd_Array1OfReal & theFlatKnots) const
788 {
789   Standard_Integer num_poles = theData.myNbPnt + 3;
790   TColStd_Array1OfInteger ContactOrder(1,num_poles);
791   TColStd_Array1OfReal    aParameters(1, num_poles);
792
793   // Fill tables taking attention to end values.
794   ContactOrder.Init(0);
795   ContactOrder(2) = ContactOrder(num_poles - 1) = 1;
796
797   theFlatKnots(1) = theFlatKnots(2) = theFlatKnots(3) = theFlatKnots(4) = theData.myC3dPF;
798   theFlatKnots(num_poles + 1) = theFlatKnots(num_poles + 2) = 
799     theFlatKnots(num_poles + 3) = theFlatKnots(num_poles + 4) = theData.myC3dPL;
800
801   thePoles(1) = theData.myC2dPF; thePoles(num_poles) = theData.myC2dPL;
802   thePoles(2) = aTangFirst; thePoles(num_poles - 1) = aTangLast;
803
804   aParameters(1) = aParameters(2) = theData.myC3dPF;
805   aParameters(num_poles - 1) = aParameters(num_poles) = theData.myC3dPL;
806
807   for (Standard_Integer ii = 3; ii <= num_poles - 2; ii++)
808   {
809     thePoles(ii) = theData.myPC2d[ii - 2];
810     aParameters(ii) = theFlatKnots(ii+2) = theData.myPC3d[ii - 2];
811   }
812   Standard_Integer inversion_problem;
813   BSplCLib::Interpolate(3,theFlatKnots,aParameters,ContactOrder,
814                         1,thePoles(1),inversion_problem);
815   if(inversion_problem)
816   {
817     return Standard_False;
818   }
819
820   return Standard_True;
821 }
822
823 //=======================================================================
824 //function : IncreaseNbPoles
825 //purpose  : Sub-method in Build.
826 //=======================================================================
827 Standard_Boolean Approx_SameParameter::IncreaseNbPoles(const TColStd_Array1OfReal & thePoles,
828                                                        const TColStd_Array1OfReal & theFlatKnots,
829                                                        Approx_SameParameter_Data & theData,
830                                                        Standard_Real &theBestSqTol) const
831 {
832   Extrema_LocateExtPC Projector;
833   Projector.Initialize(myC3d->Curve(), myC3d->FirstParameter(), myC3d->LastParameter(), theData.myTol);
834   Standard_Real curp = 0.0;
835   Standard_Boolean projok = Standard_False;
836
837   // Project middle point to fix parameterization and check projection existence.
838   const Standard_Integer aDegree = 3;
839   const Standard_Integer DerivativeRequest = 0;
840   Standard_Integer extrap_mode[2] = {aDegree, aDegree};
841   Standard_Real eval_result;
842   Standard_Real *PolesArray = (Standard_Real *) &thePoles(thePoles.Lower());
843   Standard_Integer newcount = 0;
844   for (Standard_Integer ii = 0; ii < theData.myNbPnt; ii++)
845   {
846     theData.myNewPC2d[newcount] = theData.myPC2d[ii];
847     theData.myNewPC3d[newcount] = theData.myPC3d[ii];
848     newcount++;
849
850     if(theData.myNbPnt - ii + newcount == myMaxArraySize) continue;
851
852     BSplCLib::Eval(0.5*(theData.myPC3d[ii]+theData.myPC3d[ii+1]), Standard_False, DerivativeRequest,
853       extrap_mode[0], 3, theFlatKnots, 1, PolesArray[0], eval_result);
854
855     if(eval_result < theData.myPC2d[ii] || eval_result > theData.myPC2d[ii+1])
856     {
857       Standard_Real ucons = 0.5*(theData.myPC2d[ii]+theData.myPC2d[ii+1]);
858       Standard_Real uc3d  = 0.5*(theData.myPC3d[ii]+theData.myPC3d[ii+1]);
859
860       gp_Pnt Pcons;
861       theData.myCOnS.D0(ucons,Pcons);
862       Projector.Perform(Pcons, uc3d);
863       if (Projector.IsDone())
864       {
865         curp = Projector.Point().Parameter();
866         Standard_Real dist_2 = Projector.SquareDistance();
867         if(dist_2 > theBestSqTol) theBestSqTol = dist_2;
868         projok = 1;
869       }
870       else
871       {
872         ProjectPointOnCurve(uc3d,Pcons,theData.myTol,30,myC3d->Curve(),projok,curp);
873       }
874       if(projok)
875       {
876         if(curp > theData.myPC3d[ii] + myDeltaMin && curp < theData.myPC3d[ii+1] - myDeltaMin)
877         {
878           theData.myNewPC3d[newcount] = curp;
879           theData.myNewPC2d[newcount] = ucons;
880           newcount ++;
881         }
882       }
883     }
884   } //     for (ii = 0; ii < count; ii++)
885   theData.myNewPC3d[newcount] = theData.myPC3d[theData.myNbPnt];
886   theData.myNewPC2d[newcount] = theData.myPC2d[theData.myNbPnt];
887
888   if((theData.myNbPnt != newcount) && newcount < myMaxArraySize - 1)
889   {
890     // Distribution is changed.
891     theData.Swap(newcount);
892     return Standard_True;
893   }
894
895   // Increase number of samples in two times.
896   newcount = 0;
897   for(Standard_Integer n = 0; n < theData.myNbPnt; n++)
898   {
899     theData.myNewPC3d[newcount] = theData.myPC3d[n];
900     theData.myNewPC2d[newcount] = theData.myPC2d[n];
901     newcount ++;
902
903     if(theData.myNbPnt - n + newcount == myMaxArraySize) continue;
904
905     Standard_Real ucons = 0.5*(theData.myPC2d[n]+theData.myPC2d[n+1]);
906     Standard_Real uc3d  = 0.5*(theData.myPC3d[n]+theData.myPC3d[n+1]);
907
908     gp_Pnt Pcons;
909     theData.myCOnS.D0(ucons,Pcons);
910     Projector.Perform(Pcons, uc3d);
911     if (Projector.IsDone()) 
912     {
913       curp = Projector.Point().Parameter();
914       Standard_Real dist_2 = Projector.SquareDistance();
915       if(dist_2 > theBestSqTol) theBestSqTol = dist_2;
916       projok = 1;
917     }
918     else 
919     {
920       ProjectPointOnCurve(uc3d,Pcons,theData.myTol,30,myC3d->Curve(),projok,curp);
921     }
922     if(projok)
923     {
924       if(curp > theData.myPC3d[n] + myDeltaMin && curp < theData.myPC3d[n+1] - myDeltaMin)
925       {
926         theData.myNewPC3d[newcount] = curp;
927         theData.myNewPC2d[newcount] = ucons;
928         newcount ++;
929       }
930     }
931   }
932   theData.myNewPC3d[newcount] = theData.myPC3d[theData.myNbPnt];
933   theData.myNewPC2d[newcount] = theData.myPC2d[theData.myNbPnt];
934
935   if(theData.myNbPnt != newcount)
936   {
937     // Distribution is changed.
938     theData.Swap(newcount);
939     return Standard_True;
940   }
941
942   return Standard_False;
943 }