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 }