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