src/GeomAPI/GeomAPI_Interpolate.cxx

   1 // Created on: 1994-08-18
   2 // Created by: Laurent PAINNOT
   3 // Copyright (c) 1994-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 // 08-Aug-95 : xab : interpolation uses BSplCLib::Interpolate
  18
  19 #include <BSplCLib.hxx>
  20 #include <GeomAPI_Interpolate.hxx>
  21 #include <gp_Vec.hxx>
  22 #include <PLib.hxx>
  23 #include <Standard_ConstructionError.hxx>
  24 #include <StdFail_NotDone.hxx>
  25 #include <TColgp_Array1OfPnt.hxx>
  26 #include <TColgp_Array1OfVec.hxx>
  27 #include <TColStd_Array1OfBoolean.hxx>
  28 #include <TColStd_Array1OfInteger.hxx>
  29 #include <TColStd_HArray1OfReal.hxx>
  30
  31 //=======================================================================
  32 //function : CheckPoints
  33 //purpose  :
  34 //=======================================================================
  35 static Standard_Boolean CheckPoints(const TColgp_Array1OfPnt& PointArray,
  36                                     const Standard_Real    Tolerance)
  37 {
  38   Standard_Integer ii ;
  39   Standard_Real tolerance_squared = Tolerance * Tolerance,
  40   distance_squared ;
  41   Standard_Boolean result = Standard_True ;
  42   for (ii = PointArray.Lower() ; result && ii < PointArray.Upper() ; ii++) {
  43     distance_squared =
  44       PointArray.Value(ii).SquareDistance(PointArray.Value(ii+1)) ;
  45     result = (distance_squared >= tolerance_squared) ;
  46   }
  47  return result ;
  48
  49 }
  50 //=======================================================================
  51 //function : CheckTangents
  52 //purpose  :
  53 //=======================================================================
  54 static Standard_Boolean CheckTangents(
  55                                 const TColgp_Array1OfVec&      Tangents,
  56                                 const TColStd_Array1OfBoolean& TangentFlags,
  57                                 const Standard_Real            Tolerance)
  58 {
  59   Standard_Integer ii,
  60   index ;
  61   Standard_Real tolerance_squared = Tolerance * Tolerance,
  62   distance_squared ;
  63   Standard_Boolean result = Standard_True ;
  64   index = TangentFlags.Lower() ;
  65   for (ii = Tangents.Lower(); result && ii <= Tangents.Upper() ; ii++) {
  66     if(TangentFlags.Value(index)) {
  67       distance_squared =
  68         Tangents.Value(ii).SquareMagnitude() ;
  69       result = (distance_squared >= tolerance_squared) ;
  70     }
  71     index += 1 ;
  72   }
  73  return result ;
  74
  75 }
  76 //=======================================================================
  77 //function : CheckParameters
  78 //purpose  :
  79 //=======================================================================
  80 static Standard_Boolean CheckParameters(const
  81                                         TColStd_Array1OfReal&   Parameters)
  82 {
  83   Standard_Integer ii ;
  84   Standard_Real distance ;
  85   Standard_Boolean result = Standard_True ;
  86   for (ii = Parameters.Lower() ; result && ii < Parameters.Upper() ; ii++) {
  87     distance =
  88       Parameters.Value(ii+1) - Parameters.Value(ii) ;
  89     result = (distance >= RealSmall()) ;
  90   }
  91  return result ;
  92 }
  93 //=======================================================================
  94 //function : BuildParameters
  95 //purpose  :
  96 //=======================================================================
  97 static void  BuildParameters(const Standard_Boolean        PeriodicFlag,
  98                              const TColgp_Array1OfPnt&     PointsArray,
  99                              Handle(TColStd_HArray1OfReal)& ParametersPtr)
 100 {
 101   Standard_Integer ii,
 102   index ;
 103   Standard_Real distance ;
 104   Standard_Integer
 105     num_parameters = PointsArray.Length() ;
 106   if (PeriodicFlag) {
 107     num_parameters += 1 ;
 108   }
 109   ParametersPtr =
 110     new TColStd_HArray1OfReal(1,
 111                               num_parameters) ;
 112   ParametersPtr->SetValue(1,0.0e0) ;
 113   index = 2 ;
 114   for (ii = PointsArray.Lower() ; ii < PointsArray.Upper() ; ii++) {
 115     distance =
 116       PointsArray.Value(ii).Distance(PointsArray.Value(ii+1)) ;
 117     ParametersPtr->SetValue(index,
 118                             ParametersPtr->Value(ii) + distance) ;
 119     index += 1 ;
 120   }
 121   if (PeriodicFlag) {
 122     distance =
 123       PointsArray.Value(PointsArray.Upper()).
 124         Distance(PointsArray.Value(PointsArray.Lower())) ;
 125     ParametersPtr->SetValue(index,
 126                             ParametersPtr->Value(ii) + distance) ;
 127   }
 128 }
 129 //=======================================================================
 130 //function : BuildPeriodicTangents
 131 //purpose  :
 132 //=======================================================================
 133
 134 static void BuildPeriodicTangent(
 135                       const TColgp_Array1OfPnt&      PointsArray,
 136                       TColgp_Array1OfVec&            TangentsArray,
 137                       TColStd_Array1OfBoolean&       TangentFlags,
 138                       const TColStd_Array1OfReal&    ParametersArray)
 139 {
 140   Standard_Integer
 141     ii,
 142     degree ;
 143   Standard_Real *point_array,
 144   *parameter_array,
 145   eval_result[2][3] ;
 146
 147   gp_Vec a_vector ;
 148
 149   if (PointsArray.Length() < 3) {
 150     throw Standard_ConstructionError();
 151     }
 152
 153   if (!TangentFlags.Value(1)) {
 154     degree = 3 ;
 155     if (PointsArray.Length() == 3) {
 156       degree = 2 ;
 157     }
 158     point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()) ;
 159     parameter_array =
 160       (Standard_Real *) &ParametersArray.Value(1) ;
 161     TangentFlags.SetValue(1,Standard_True) ;
 162     PLib::EvalLagrange(ParametersArray.Value(1),
 163                        1,
 164                        degree,
 165                        3,
 166                        point_array[0],
 167                        parameter_array[0],
 168                        eval_result[0][0]) ;
 169     for (ii = 1 ; ii <= 3 ; ii++) {
 170       a_vector.SetCoord(ii,eval_result[1][ii-1]) ;
 171     }
 172     TangentsArray.SetValue(1,a_vector) ;
 173   }
 174  }
 175 //=======================================================================
 176 //function : BuildTangents
 177 //purpose  :
 178 //=======================================================================
 179
 180 static void BuildTangents(const TColgp_Array1OfPnt&      PointsArray,
 181                           TColgp_Array1OfVec&            TangentsArray,
 182                           TColStd_Array1OfBoolean&       TangentFlags,
 183                           const TColStd_Array1OfReal&    ParametersArray)
 184 {
 185  Standard_Integer ii,
 186  degree ;
 187  Standard_Real *point_array,
 188  *parameter_array,
 189
 190  eval_result[2][3] ;
 191  gp_Vec a_vector ;
 192
 193  degree = 3 ;
 194
 195  if ( PointsArray.Length() < 3) {
 196    throw Standard_ConstructionError();
 197    }
 198  if (PointsArray.Length() == 3) {
 199    degree = 2 ;
 200  }
 201  if (!TangentFlags.Value(1)) {
 202    point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()) ;
 203    parameter_array =
 204      (Standard_Real *) &ParametersArray.Value(1) ;
 205    TangentFlags.SetValue(1,Standard_True) ;
 206    PLib::EvalLagrange(ParametersArray.Value(1),
 207                       1,
 208                       degree,
 209                       3,
 210                       point_array[0],
 211                       parameter_array[0],
 212                       eval_result[0][0]) ;
 213    for (ii = 1 ; ii <= 3 ; ii++) {
 214      a_vector.SetCoord(ii,eval_result[1][ii-1]) ;
 215    }
 216    TangentsArray.SetValue(1,a_vector) ;
 217  }
 218  if (! TangentFlags.Value(TangentFlags.Upper())) {
 219    point_array =
 220      (Standard_Real *) &PointsArray.Value(PointsArray.Upper() - degree) ;
 221    TangentFlags.SetValue(TangentFlags.Upper(),Standard_True) ;
 222    parameter_array =
 223      (Standard_Real *) &ParametersArray.Value(ParametersArray.Upper() - degree) ;
 224    PLib::EvalLagrange(ParametersArray.Value(ParametersArray.Upper()),
 225                       1,
 226                       degree,
 227                       3,
 228                       point_array[0],
 229                       parameter_array[0],
 230                       eval_result[0][0]) ;
 231    for (ii = 1 ; ii <= 3 ; ii++) {
 232      a_vector.SetCoord(ii,eval_result[1][ii-1]) ;
 233    }
 234    TangentsArray.SetValue(TangentsArray.Upper(),a_vector) ;
 235  }
 236 }
 237 //=======================================================================
 238 //function : BuildTangents
 239 //purpose  : scale the given tangent so that they have the length of
 240 // the size of the derivative of the lagrange interpolation
 241 //
 242 //=======================================================================
 243 static void ScaleTangents(const TColgp_Array1OfPnt&      PointsArray,
 244                           TColgp_Array1OfVec&            TangentsArray,
 245                           const TColStd_Array1OfBoolean& TangentFlags,
 246                           const TColStd_Array1OfReal&     ParametersArray)
 247 {
 248  Standard_Integer ii,
 249  jj,
 250  degree=0,
 251  index,
 252  num_points ;
 253
 254  Standard_Real *point_array,
 255  *parameter_array,
 256  value[2],
 257  ratio,
 258  eval_result[2][3] ;
 259
 260  gp_Vec a_vector ;
 261
 262  num_points = PointsArray.Length() ;
 263  if (num_points == 2) {
 264     degree = 1 ;
 265   }
 266  else if (num_points >= 3) {
 267    degree = 2 ;
 268  }
 269
 270  index = PointsArray.Lower() ;
 271  for (ii = TangentFlags.Lower()  ; ii <=  TangentFlags.Upper() ; ii++) {
 272    if (TangentFlags.Value(ii)) {
 273      point_array =
 274        (Standard_Real *) &PointsArray.Value(index) ;
 275      parameter_array =
 276        (Standard_Real *) &ParametersArray.Value(index) ;
 277      PLib::EvalLagrange(ParametersArray.Value(ii),
 278                         1,
 279                         degree,
 280                         3,
 281                         point_array[0],
 282                         parameter_array[0],
 283                         eval_result[0][0]) ;
 284      value[0] =
 285        value[1] = 0.0e0 ;
 286      for (jj = 1 ; jj <= 3 ; jj++) {
 287        value[0] += Abs(TangentsArray.Value(ii).Coord(jj)) ;
 288        value[1] += Abs(eval_result[1][jj-1]) ;
 289      }
 290      ratio = value[1] / value[0] ;
 291      for (jj = 1 ; jj <= 3 ; jj++) {
 292        a_vector.SetCoord(jj, ratio *
 293                          TangentsArray.Value(ii).Coord(jj)) ;
 294      }
 295      TangentsArray.SetValue(ii, a_vector) ;
 296      if (ii != TangentFlags.Lower()) {
 297        index += 1 ;
 298      }
 299      if (index >  PointsArray.Upper() - degree) {
 300        index = PointsArray.Upper() - degree ;
 301      }
 302
 303    }
 304  }
 305 }
 306
 307 //=======================================================================
 308 //function : GeomAPI_Interpolate
 309 //purpose  :
 310 //=======================================================================
 311
 312 GeomAPI_Interpolate::GeomAPI_Interpolate
 313    (const Handle(TColgp_HArray1OfPnt)& PointsPtr,
 314     const Standard_Boolean            PeriodicFlag,
 315     const Standard_Real               Tolerance) :
 316 myTolerance(Tolerance),
 317 myPoints(PointsPtr),
 318 myIsDone(Standard_False),
 319 myPeriodic(PeriodicFlag),
 320 myTangentRequest(Standard_False)
 321 {
 322  Standard_Integer ii ;
 323  Standard_Boolean result =
 324    CheckPoints(PointsPtr->Array1(),
 325                Tolerance) ;
 326  myTangents =
 327      new TColgp_HArray1OfVec(myPoints->Lower(),
 328                               myPoints->Upper()) ;
 329  myTangentFlags =
 330       new TColStd_HArray1OfBoolean(myPoints->Lower(),
 331                                    myPoints->Upper()) ;
 332
 333  if (!result) {
 334    throw Standard_ConstructionError();
 335    }
 336  BuildParameters(PeriodicFlag,
 337                  PointsPtr->Array1(),
 338                  myParameters) ;
 339
 340   for (ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) {
 341     myTangentFlags->SetValue(ii,Standard_False) ;
 342   }
 343
 344
 345
 346 }
 347
 348 //=======================================================================
 349 //function : GeomAPI_Interpolate
 350 //purpose  :
 351 //=======================================================================
 352
 353 GeomAPI_Interpolate::GeomAPI_Interpolate
 354    (const Handle(TColgp_HArray1OfPnt)&    PointsPtr,
 355     const Handle(TColStd_HArray1OfReal)&  ParametersPtr,
 356     const Standard_Boolean               PeriodicFlag,
 357     const Standard_Real                  Tolerance) :
 358 myTolerance(Tolerance),
 359 myPoints(PointsPtr),
 360 myIsDone(Standard_False),
 361 myParameters(ParametersPtr),
 362 myPeriodic(PeriodicFlag),
 363 myTangentRequest(Standard_False)
 364 {
 365  Standard_Integer ii ;
 366
 367
 368  Standard_Boolean result =
 369    CheckPoints(PointsPtr->Array1(),
 370                Tolerance) ;
 371
 372  if (PeriodicFlag) {
 373    if ((PointsPtr->Length()) + 1 != ParametersPtr->Length()) {
 374      throw Standard_ConstructionError();
 375    }
 376  }
 377  myTangents =
 378      new TColgp_HArray1OfVec(myPoints->Lower(),
 379                               myPoints->Upper()) ;
 380  myTangentFlags =
 381       new TColStd_HArray1OfBoolean(myPoints->Lower(),
 382                                     myPoints->Upper()) ;
 383
 384  if (!result) {
 385    throw Standard_ConstructionError();
 386    }
 387
 388  result =
 389  CheckParameters(ParametersPtr->Array1()) ;
 390  if (!result) {
 391    throw Standard_ConstructionError();
 392    }
 393
 394  for (ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) {
 395    myTangentFlags->SetValue(ii,Standard_False) ;
 396  }
 397
 398 }
 399 //=======================================================================
 400 //function : Load
 401 //purpose  :
 402 //=======================================================================
 403
 404 void GeomAPI_Interpolate::Load(
 405    const TColgp_Array1OfVec&              Tangents,
 406    const Handle(TColStd_HArray1OfBoolean)& TangentFlagsPtr,
 407    const Standard_Boolean Scale )
 408
 409 {
 410  Standard_Boolean result ;
 411  Standard_Integer ii ;
 412  myTangentRequest = Standard_True ;
 413  myTangentFlags = TangentFlagsPtr ;
 414  if (Tangents.Length() != myPoints->Length() ||
 415      TangentFlagsPtr->Length() != myPoints->Length()) {
 416    throw Standard_ConstructionError();
 417    }
 418  result  =
 419    CheckTangents(Tangents,
 420                  TangentFlagsPtr->Array1(),
 421                  myTolerance) ;
 422  if (result) {
 423     myTangents =
 424       new TColgp_HArray1OfVec(Tangents.Lower(),Tangents.Upper()) ;
 425     for (ii = Tangents.Lower() ; ii <= Tangents.Upper() ; ii++ ) {
 426       myTangents->SetValue(ii,Tangents.Value(ii)) ;
 427     }
 428
 429     if(Scale) {
 430       ScaleTangents(myPoints->Array1(),
 431                     myTangents->ChangeArray1(),
 432                     TangentFlagsPtr->Array1(),
 433                     myParameters->Array1()) ;
 434     }
 435   }
 436  else {
 437    throw Standard_ConstructionError();
 438    }
 439
 440
 441 }
 442
 443 //=======================================================================
 444 //function : Load
 445 //purpose  :
 446 //=======================================================================
 447
 448 void GeomAPI_Interpolate::Load(const gp_Vec& InitialTangent,
 449                                const gp_Vec& FinalTangent,
 450                                const Standard_Boolean Scale )
 451 {
 452   Standard_Boolean result ;
 453   myTangentRequest = Standard_True ;
 454   myTangentFlags->SetValue(1,Standard_True) ;
 455   myTangentFlags->SetValue(myPoints->Length(),Standard_True) ;
 456   myTangents->SetValue(1,InitialTangent) ;
 457   myTangents->SetValue(myPoints->Length(),FinalTangent) ;
 458   result =
 459     CheckTangents(myTangents->Array1(),
 460                   myTangentFlags->Array1(),
 461                   myTolerance) ;
 462   if (!result) {
 463     throw Standard_ConstructionError();
 464     }
 465
 466   if(Scale) {
 467     ScaleTangents(myPoints->Array1(),
 468                   myTangents->ChangeArray1(),
 469                   myTangentFlags->Array1(),
 470                   myParameters->Array1()) ;
 471   }
 472
 473 }
 474 //=======================================================================
 475 //function : Perform
 476 //purpose  :
 477 //=======================================================================
 478
 479 void GeomAPI_Interpolate::Perform()
 480 {
 481   if (myPeriodic) {
 482     PerformPeriodic() ;
 483   }
 484   else {
 485     PerformNonPeriodic() ;
 486   }
 487 }
 488 //=======================================================================
 489 //function : PerformPeriodic
 490 //purpose  :
 491 //=======================================================================
 492
 493 void GeomAPI_Interpolate::PerformPeriodic()
 494 {
 495   Standard_Integer degree,
 496   ii,
 497   jj,
 498   index,
 499   index1,
 500 //  index2,
 501   mult_index,
 502   half_order,
 503   inversion_problem,
 504   num_points,
 505   num_distinct_knots,
 506   num_poles  ;
 507
 508   Standard_Real period ;
 509
 510   gp_Pnt a_point ;
 511
 512   num_points = myPoints->Length() ;
 513   period = myParameters->Value(myParameters->Upper()) -
 514    myParameters->Value(myParameters->Lower())  ;
 515   num_poles = num_points + 1 ;
 516   if (num_points == 2 && !myTangentRequest) {
 517 //
 518 // build a periodic curve of degree 1
 519 //
 520
 521     degree = 1 ;
 522     TColStd_Array1OfInteger  deg1_mults(1,num_poles) ;
 523     for (ii = 1 ; ii <= num_poles ; ii++) {
 524       deg1_mults.SetValue(ii,1) ;
 525     }
 526
 527     myCurve =
 528       new Geom_BSplineCurve(myPoints->Array1(),
 529                             myParameters->Array1(),
 530                             deg1_mults,
 531                             degree,
 532                             myPeriodic) ;
 533     myIsDone = Standard_True ;
 534
 535   }
 536   else {
 537     num_distinct_knots = num_points + 1 ;
 538     half_order = 2 ;
 539     degree = 3 ;
 540     num_poles += 2 ;
 541     if (myTangentRequest)
 542       for (ii = myTangentFlags->Lower() + 1 ;
 543            ii <= myTangentFlags->Upper() ; ii++) {
 544         if (myTangentFlags->Value(ii)) {
 545           num_poles += 1 ;
 546         }
 547       }
 548
 549     TColStd_Array1OfReal     parameters(1,num_poles) ;
 550     TColStd_Array1OfReal     flatknots(1,num_poles + degree + 1) ;
 551     TColStd_Array1OfInteger  mults(1,num_distinct_knots) ;
 552     TColStd_Array1OfInteger  contact_order_array(1, num_poles) ;
 553     TColgp_Array1OfPnt       poles(1,num_poles) ;
 554
 555     for (ii = 1 ; ii <= half_order ; ii++) {
 556       flatknots.SetValue(ii,myParameters->Value(myParameters->Upper() -1) -
 557                          period) ;
 558       flatknots.SetValue(ii + half_order,myParameters->
 559                          Value(myParameters->Lower())) ;
 560       flatknots.SetValue(num_poles + ii,
 561                          myParameters->Value(myParameters->Upper())) ;
 562       flatknots.SetValue(num_poles + half_order + ii,
 563                          myParameters->Value(half_order) + period)  ;
 564     }
 565     for (ii = 1 ; ii <= num_poles ; ii++) {
 566       contact_order_array.SetValue(ii,0)  ;
 567     }
 568     for (ii = 2 ; ii < num_distinct_knots ; ii++) {
 569       mults.SetValue(ii,1) ;
 570     }
 571     mults.SetValue(1,half_order) ;
 572     mults.SetValue(num_distinct_knots ,half_order) ;
 573     if (num_points >= 3) {
 574
 575 //
 576 //   only enter here if there are more than 3 points otherwise
 577 //   it means we have already the tangent
 578 //
 579       BuildPeriodicTangent(myPoints->Array1(),
 580                            myTangents->ChangeArray1(),
 581                            myTangentFlags->ChangeArray1(),
 582                            myParameters->Array1()) ;
 583     }
 584     contact_order_array.SetValue(2,1)  ;
 585     parameters.SetValue(1,myParameters->Value(1)) ;
 586     parameters.SetValue(2,myParameters->Value(1)) ;
 587     poles.SetValue(1,myPoints->Value(1)) ;
 588     for (jj = 1 ; jj <= 3 ; jj++) {
 589       a_point.SetCoord(jj,myTangents->Value(1).Coord(jj)) ;
 590     }
 591     poles.SetValue(2,a_point) ;
 592
 593     mult_index = 2 ;
 594     index = 3 ;
 595     index1 = degree + 2 ;
 596     if (myTangentRequest) {
 597       for (ii = myTangentFlags->Lower() + 1 ;
 598            ii <= myTangentFlags->Upper() ; ii++) {
 599         parameters.SetValue(index,myParameters->Value(ii)) ;
 600         flatknots.SetValue(index1,myParameters->Value(ii)) ;
 601         poles.SetValue(index,myPoints->Value(ii)) ;
 602         index += 1  ;
 603         index1 += 1 ;
 604         if (myTangentFlags->Value(ii)) {
 605           mults.SetValue(mult_index,mults.Value(mult_index)  + 1)   ;
 606           contact_order_array(index) = 1 ;
 607
 608           parameters.SetValue(index,
 609                               myParameters->Value(ii)) ;
 610           flatknots.SetValue(index1,myParameters->Value(ii)) ;
 611           for (jj = 1 ; jj <= 3 ; jj++) {
 612             a_point.SetCoord(jj,myTangents->Value(ii).Coord(jj)) ;
 613           }
 614           poles.SetValue(index,a_point) ;
 615           index  += 1 ;
 616           index1 += 1 ;
 617         }
 618        mult_index += 1 ;
 619       }
 620     }
 621     else {
 622       index = degree + 1 ;
 623       index1 = 2 ;
 624       for(ii = myParameters->Lower()  ; ii <= myParameters->Upper()  ; ii++) {
 625         parameters.SetValue(index1,
 626                             myParameters->Value(ii)) ;
 627         flatknots.SetValue(index,
 628                            myParameters->Value(ii)) ;
 629         index += 1 ;
 630         index1 += 1 ;
 631       }
 632       index = 3 ;
 633       for (ii = myPoints->Lower() + 1 ; ii <= myPoints->Upper() ; ii++) {
 634 //
 635 // copy all the given points since the last one will be initialized
 636 // below by the first point in the array myPoints
 637 //
 638         poles.SetValue(index,
 639                        myPoints->Value(ii)) ;
 640         index += 1 ;
 641       }
 642
 643     }
 644     contact_order_array.SetValue(num_poles - 1, 1)  ;
 645     parameters.SetValue(num_poles-1,
 646                         myParameters->Value(myParameters->Upper())) ;
 647 //
 648 // for the periodic curve ONLY  the tangent of the first point
 649 // will be used since the curve should close itself at the first
 650 // point See BuildPeriodicTangent
 651 //
 652     for (jj = 1 ; jj <= 3 ; jj++) {
 653       a_point.SetCoord(jj,myTangents->Value(1).Coord(jj)) ;
 654     }
 655     poles.SetValue(num_poles-1,a_point) ;
 656
 657     parameters.SetValue(num_poles,
 658                         myParameters->Value(myParameters->Upper())) ;
 659
 660     poles.SetValue(num_poles,
 661                      myPoints->Value(1)) ;
 662
 663
 664     BSplCLib::Interpolate(degree,
 665                           flatknots,
 666                           parameters,
 667                           contact_order_array,
 668                           poles,
 669                           inversion_problem) ;
 670     if (!inversion_problem) {
 671       TColgp_Array1OfPnt   newpoles(poles.Value(1),
 672                                     1,
 673                                     num_poles - 2) ;
 674       myCurve =
 675         new Geom_BSplineCurve(newpoles,
 676                               myParameters->Array1(),
 677                               mults,
 678                               degree,
 679                               myPeriodic) ;
 680       myIsDone = Standard_True ;
 681   }
 682  }
 683 }
 684
 685
 686 //=======================================================================
 687 //function : PerformNonPeriodic
 688 //purpose  :
 689 //=======================================================================
 690
 691 void GeomAPI_Interpolate::PerformNonPeriodic()
 692 {
 693   Standard_Integer degree,
 694   ii,
 695   jj,
 696   index,
 697   index1,
 698   index2,
 699   index3,
 700   mult_index,
 701   inversion_problem,
 702   num_points,
 703   num_distinct_knots,
 704   num_poles  ;
 705
 706   gp_Pnt a_point ;
 707
 708   num_points =
 709   num_distinct_knots =
 710   num_poles = myPoints->Length() ;
 711   if (num_poles == 2 &&   !myTangentRequest)  {
 712     degree = 1 ;
 713   }
 714   else if (num_poles == 3 && !myTangentRequest) {
 715     degree = 2 ;
 716     num_distinct_knots = 2 ;
 717   }
 718   else {
 719     degree = 3 ;
 720     num_poles += 2 ;
 721     if (myTangentRequest)
 722       for (ii = myTangentFlags->Lower() + 1 ;
 723            ii < myTangentFlags->Upper() ; ii++) {
 724         if (myTangentFlags->Value(ii)) {
 725           num_poles += 1 ;
 726         }
 727       }
 728     }
 729
 730
 731   TColStd_Array1OfReal     parameters(1,num_poles) ;
 732   TColStd_Array1OfReal     flatknots(1,num_poles + degree + 1) ;
 733   TColStd_Array1OfInteger  mults(1,num_distinct_knots) ;
 734   TColStd_Array1OfReal     knots(1,num_distinct_knots) ;
 735   TColStd_Array1OfInteger  contact_order_array(1, num_poles) ;
 736   TColgp_Array1OfPnt       poles(1,num_poles) ;
 737
 738   for (ii = 1 ; ii <= degree + 1 ; ii++) {
 739     flatknots.SetValue(ii,myParameters->Value(1)) ;
 740     flatknots.SetValue(ii + num_poles,
 741                        myParameters->Value(num_points)) ;
 742   }
 743   for (ii = 1 ; ii <= num_poles ; ii++) {
 744     contact_order_array.SetValue(ii,0)  ;
 745   }
 746   for (ii = 2 ; ii < num_distinct_knots ; ii++) {
 747     mults.SetValue(ii,1) ;
 748   }
 749   mults.SetValue(1,degree + 1) ;
 750   mults.SetValue(num_distinct_knots ,degree + 1) ;
 751
 752   switch (degree) {
 753   case 1:
 754     for (ii = 1 ; ii <= num_poles ; ii++) {
 755       poles.SetValue(ii ,myPoints->Value(ii)) ;
 756     }
 757     myCurve =
 758       new Geom_BSplineCurve(poles,
 759                             myParameters->Array1(),
 760                             mults,
 761                             degree) ;
 762     myIsDone = Standard_True ;
 763     break ;
 764   case 2:
 765     knots.SetValue(1,myParameters->Value(1)) ;
 766     knots.SetValue(2,myParameters->Value(3)) ;
 767     for (ii = 1 ; ii <= num_poles ; ii++) {
 768       poles.SetValue(ii,myPoints->Value(ii)) ;
 769
 770     }
 771     BSplCLib::Interpolate(degree,
 772                           flatknots,
 773                           myParameters->Array1(),
 774                           contact_order_array,
 775                           poles,
 776                           inversion_problem) ;
 777     if (!inversion_problem) {
 778       myCurve =
 779         new Geom_BSplineCurve(poles,
 780                               knots,
 781                               mults,
 782                               degree) ;
 783       myIsDone = Standard_True ;
 784     }
 785     break ;
 786   case 3:
 787 //
 788 // check if the boundary conditions are set
 789 //
 790     if (num_points >= 3) {
 791 //
 792 // cannot build the tangents with degree 3 with only 2 points
 793 // if those where not given in advance
 794 //
 795       BuildTangents(myPoints->Array1(),
 796                     myTangents->ChangeArray1(),
 797                     myTangentFlags->ChangeArray1(),
 798                     myParameters->Array1()) ;
 799     }
 800     contact_order_array.SetValue(2,1)  ;
 801     parameters.SetValue(1,myParameters->Value(1)) ;
 802     parameters.SetValue(2,myParameters->Value(1)) ;
 803     poles.SetValue(1,myPoints->Value(1)) ;
 804     for (jj = 1 ; jj <= 3 ; jj++) {
 805       a_point.SetCoord(jj,myTangents->Value(1).Coord(jj)) ;
 806
 807     }
 808     poles.SetValue(2,a_point) ;
 809     mult_index = 2 ;
 810     index = 3 ;
 811     index1 = 2 ;
 812     index2 = myPoints->Lower() + 1 ;
 813     index3 = degree + 2 ;
 814     if (myTangentRequest) {
 815       for (ii = myParameters->Lower() + 1 ;
 816            ii < myParameters->Upper() ; ii++) {
 817         parameters.SetValue(index,myParameters->Value(ii)) ;
 818         poles.SetValue(index,myPoints->Value(index2)) ;
 819         flatknots.SetValue(index3,myParameters->Value(ii)) ;
 820         index += 1 ;
 821         index3 += 1 ;
 822         if (myTangentFlags->Value(index1)) {
 823 //
 824 // set the multiplicities, the order of the contact, the
 825 // the flatknots,
 826 //
 827           mults.SetValue(mult_index,mults.Value(mult_index) + 1) ;
 828           contact_order_array(index) = 1 ;
 829           flatknots.SetValue(index3, myParameters->Value(ii)) ;
 830           parameters.SetValue(index,
 831                               myParameters->Value(ii)) ;
 832           for (jj = 1 ; jj <= 3 ; jj++) {
 833             a_point.SetCoord(jj,myTangents->Value(ii).Coord(jj)) ;
 834           }
 835           poles.SetValue(index,a_point) ;
 836           index += 1  ;
 837           index3 += 1 ;
 838         }
 839         mult_index += 1 ;
 840         index1 += 1 ;
 841         index2 += 1 ;
 842
 843       }
 844     }
 845     else {
 846       index1 = 2 ;
 847       for(ii = myParameters->Lower()  ; ii <= myParameters->Upper()  ; ii++) {
 848         parameters.SetValue(index1,
 849                             myParameters->Value(ii)) ;
 850         index1 += 1 ;
 851       }
 852       index = 3 ;
 853       for (ii = myPoints->Lower() + 1 ; ii <= myPoints->Upper() - 1 ; ii++) {
 854         poles.SetValue(index,
 855                        myPoints->Value(ii)) ;
 856         index += 1 ;
 857       }
 858
 859
 860       index = degree + 1 ;
 861       for(ii = myParameters->Lower()  ; ii <= myParameters->Upper()  ; ii++) {
 862         flatknots.SetValue(index,
 863                            myParameters->Value(ii)) ;
 864         index += 1 ;
 865       }
 866     }
 867     for (jj = 1 ; jj <= 3 ; jj++) {
 868       a_point.SetCoord(jj,
 869                        myTangents->Value(num_points).Coord(jj)) ;
 870     }
 871     poles.SetValue(num_poles-1 ,a_point) ;
 872
 873     contact_order_array.SetValue(num_poles - 1,1) ;
 874     parameters.SetValue(num_poles,
 875                         myParameters->Value(myParameters->Upper())) ;
 876     parameters.SetValue(num_poles -1,
 877                         myParameters->Value(myParameters->Upper())) ;
 878
 879     poles.SetValue(num_poles,
 880                    myPoints->Value(num_points)) ;
 881
 882     BSplCLib::Interpolate(degree,
 883                           flatknots,
 884                           parameters,
 885                           contact_order_array,
 886                           poles,
 887                           inversion_problem) ;
 888     if (!inversion_problem) {
 889       myCurve =
 890         new Geom_BSplineCurve(poles,
 891                               myParameters->Array1(),
 892                               mults,
 893                               degree) ;
 894       myIsDone = Standard_True ;
 895     }
 896     break ;
 897
 898   }
 899 }
 900 //=======================================================================
 901 //function : Handle(Geom_BSplineCurve)&
 902 //purpose  :
 903 //=======================================================================
 904
 905 const Handle(Geom_BSplineCurve)& GeomAPI_Interpolate::Curve() const
 906 {
 907   if ( !myIsDone)
 908     throw StdFail_NotDone(" ");
 909   return myCurve;
 910 }
 911
 912
 913
 914 //=======================================================================
 915 //function : Geom_BSplineCurve
 916 //purpose  :
 917 //=======================================================================
 918
 919 GeomAPI_Interpolate::operator Handle(Geom_BSplineCurve)() const
 920 {
 921   return myCurve;
 922 }
 923
 924
 925 //=======================================================================
 926 //function : IsDone
 927 //purpose  :
 928 //=======================================================================
 929
 930 Standard_Boolean GeomAPI_Interpolate::IsDone() const
 931 {
 932   return myIsDone;
 933 }
 934