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