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