src/GeomAPI/GeomAPI_PointsToBSplineSurface.cxx

   1 // Created on: 1995-01-17
   2 // Created by: Remi LEQUETTE
   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 <AppDef_BSplineCompute.hxx>
  19 #include <AppDef_Variational.hxx>
  20 #include <AppParCurves_ConstraintCouple.hxx>
  21 #include <AppParCurves_HArray1OfConstraintCouple.hxx>
  22 #include <BSplCLib.hxx>
  23 #include <Geom_BSplineCurve.hxx>
  24 #include <Geom_BSplineSurface.hxx>
  25 #include <GeomAPI_PointsToBSpline.hxx>
  26 #include <GeomAPI_PointsToBSplineSurface.hxx>
  27 #include <GeomFill_AppSurf.hxx>
  28 #include <GeomFill_Line.hxx>
  29 #include <GeomFill_SectionGenerator.hxx>
  30 #include <gp_Pnt.hxx>
  31 #include <math_Vector.hxx>
  32 #include <StdFail_NotDone.hxx>
  33 #include <TColgp_Array1OfPnt.hxx>
  34 #include <AppDef_BSpParLeastSquareOfMyBSplGradientOfBSplineCompute.hxx>
  35
  36 static void BuildParameters(const AppDef_MultiLine& theLine,
  37   const Approx_ParametrizationType theParT,
  38   TColStd_Array1OfReal&  thePars)
  39 {
  40   Standard_Integer i, j, nbP3d = theLine.NbPoints();
  41   Standard_Real dist;
  42   Standard_Integer firstP = 1, lastP = theLine.NbMultiPoints();
  43   const Standard_Integer aNbp = lastP - firstP + 1;
  44
  45
  46   if (aNbp == 2) {
  47     thePars(firstP) = 0.0;
  48     thePars(lastP) = 1.0;
  49   }
  50   else if (theParT == Approx_ChordLength || theParT == Approx_Centripetal)
  51   {
  52
  53     thePars(firstP) = 0.0;
  54     dist = 0.0;
  55
  56     for (i = firstP + 1; i <= lastP; i++)
  57     {
  58       AppDef_MultiPointConstraint aMPC = theLine.Value(i - 1);
  59       AppDef_MultiPointConstraint aMPC1 = theLine.Value(i);
  60
  61       dist = 0.0;
  62       for (j = 1; j <= nbP3d; j++)
  63       {
  64         const gp_Pnt &aP1 = aMPC.Point(j),
  65                      &aP2 = aMPC1.Point(j);
  66         dist += aP2.SquareDistance(aP1);
  67       }
  68
  69       dist = Sqrt(dist);
  70       if (theParT == Approx_ChordLength)
  71       {
  72         thePars(i) = thePars(i - 1) + dist;
  73       }
  74       else
  75       {// Par == Approx_Centripetal
  76         thePars(i) = thePars(i - 1) + Sqrt(dist);
  77       }
  78     }
  79     for (i = firstP; i <= lastP; i++) thePars(i) /= thePars(lastP);
  80   }
  81   else {
  82     for (i = firstP; i <= lastP; i++) {
  83       thePars(i) = (Standard_Real(i) - firstP) /
  84         (Standard_Real(lastP - Standard_Real(firstP)));
  85     }
  86   }
  87
  88 }
  89
  90 static void BuildPeriodicTangent(const AppDef_MultiLine& theLine,
  91   const TColStd_Array1OfReal&  thePars,
  92   math_Vector& theTang)
  93 {
  94   Standard_Integer firstpt = 1, lastpt = theLine.NbMultiPoints();
  95   Standard_Integer nbpoints = lastpt - firstpt + 1;
  96   //
  97   if (nbpoints <= 2)
  98   {
  99     return;
 100   }
 101   //
 102   Standard_Integer i, nnpol, nnp = Min(nbpoints, 9);
 103   nnpol = nnp;
 104   Standard_Integer lastp = Min(lastpt, firstpt + nnp - 1);
 105   Standard_Real U;
 106   AppParCurves_Constraint Cons = AppParCurves_TangencyPoint;
 107   if (nnp <= 4)
 108   {
 109     Cons = AppParCurves_PassPoint;
 110   }
 111   Standard_Integer nbP = 3 * theLine.NbPoints();
 112   math_Vector V1(1, nbP), V2(1, nbP);
 113   math_Vector P1(firstpt, lastp);
 114   //
 115   for (i = firstpt; i <= lastp; i++)
 116   {
 117     P1(i) = thePars(i);
 118   }
 119
 120   AppDef_BSpParLeastSquareOfMyBSplGradientOfBSplineCompute SQ1(theLine, firstpt, lastp, Cons, Cons, nnpol);
 121   SQ1.Perform(P1);
 122   const AppParCurves_MultiCurve& C1 = SQ1.BezierValue();
 123   U = 0.0;
 124   Standard_Integer  j, nbP3d = theLine.NbPoints();
 125
 126   gp_Pnt aP;
 127   gp_Vec aV;
 128   j = 1;
 129   for (i = 1; i <= nbP3d; i++) {
 130     C1.D1(i, U, aP, aV);
 131     V1(j) = aV.X();
 132     V1(j + 1) = aV.Y();
 133     V1(j + 2) = aV.Z();
 134     j += 3;
 135   }
 136
 137   Standard_Integer firstp = Max(firstpt, lastpt - nnp + 1);
 138
 139   if (firstp == firstpt && lastp == lastpt) {
 140     U = 1.0;
 141     j = 1;
 142     for (i = 1; i <= nbP3d; i++) {
 143       C1.D1(i, U, aP, aV);
 144       V2(j) = aV.X();
 145       V2(j + 1) = aV.Y();
 146       V2(j + 2) = aV.Z();
 147       j += 3;
 148     }
 149   }
 150   else {
 151     AppDef_BSpParLeastSquareOfMyBSplGradientOfBSplineCompute
 152       SQ2(theLine, firstp, lastpt, Cons, Cons, nnpol);
 153
 154     math_Vector P2(firstp, lastpt);
 155     for (i = firstp; i <= lastpt; i++) P2(i) = thePars(i);
 156     SQ2.Perform(P2);
 157
 158     const AppParCurves_MultiCurve& C2 = SQ2.BezierValue();
 159     U = 1.0;
 160     j = 1;
 161     for (i = 1; i <= nbP3d; i++) {
 162       C2.D1(i, U, aP, aV);
 163       V2(j) = aV.X();
 164       V2(j + 1) = aV.Y();
 165       V2(j + 2) = aV.Z();
 166       j += 3;
 167     }
 168   }
 169
 170   theTang = 0.5*(V1 + V2);
 171
 172 }
 173
 174 //=======================================================================
 175 //function : GeomAPI_PointsToBSplineSurface
 176 //purpose  :
 177 //=======================================================================
 178 GeomAPI_PointsToBSplineSurface::GeomAPI_PointsToBSplineSurface()
 179 : myIsDone ( Standard_False)
 180 {
 181 }
 182
 183
 184 //=======================================================================
 185 //function : GeomAPI_PointsToBSplineSurface
 186 //purpose  :
 187 //=======================================================================
 188
 189 GeomAPI_PointsToBSplineSurface::GeomAPI_PointsToBSplineSurface
 190 (const TColgp_Array2OfPnt& Points,
 191  const Standard_Integer DegMin,
 192  const Standard_Integer DegMax,
 193  const GeomAbs_Shape Continuity,
 194  const Standard_Real Tol3D)
 195 : myIsDone ( Standard_False)
 196 {
 197   Init(Points,DegMin,DegMax,Continuity,Tol3D);
 198 }
 199
 200 //=======================================================================
 201 //function : GeomAPI_PointsToBSplineSurface
 202 //purpose  :
 203 //=======================================================================
 204
 205 GeomAPI_PointsToBSplineSurface::GeomAPI_PointsToBSplineSurface
 206 (const TColgp_Array2OfPnt& Points,
 207  const Approx_ParametrizationType ParType,
 208  const Standard_Integer DegMin,
 209  const Standard_Integer DegMax,
 210  const GeomAbs_Shape Continuity,
 211  const Standard_Real Tol3D)
 212 : myIsDone ( Standard_False)
 213 {
 214   Init(Points,ParType,DegMin,DegMax,Continuity,Tol3D);
 215 }
 216
 217 //=======================================================================
 218 //function : GeomAPI_PointsToBSplineSurface
 219 //purpose  :
 220 //=======================================================================
 221
 222 GeomAPI_PointsToBSplineSurface::GeomAPI_PointsToBSplineSurface
 223 (const TColgp_Array2OfPnt& Points,
 224  const Standard_Real Weight1,
 225  const Standard_Real Weight2,
 226  const Standard_Real Weight3,
 227  const Standard_Integer DegMax,
 228  const GeomAbs_Shape Continuity,
 229  const Standard_Real Tol3D)
 230 : myIsDone ( Standard_False)
 231 {
 232   Init(Points,Weight1,Weight2,Weight3,DegMax,Continuity,Tol3D);
 233 }
 234
 235
 236
 237 //=======================================================================
 238 //function : GeomAPI_PointsToBSplineSurface
 239 //purpose  :
 240 //=======================================================================
 241
 242 GeomAPI_PointsToBSplineSurface::GeomAPI_PointsToBSplineSurface
 243 (const TColStd_Array2OfReal& ZPoints,
 244  const Standard_Real X0,
 245  const Standard_Real dX,
 246  const Standard_Real Y0,
 247  const Standard_Real dY,
 248  const Standard_Integer DegMin,
 249  const Standard_Integer DegMax,
 250  const GeomAbs_Shape Continuity,
 251  const Standard_Real Tol3D)
 252 : myIsDone ( Standard_False)
 253 {
 254   Init(ZPoints,X0,dX,Y0,dY,DegMin,DegMax,Continuity,Tol3D);
 255 }
 256
 257
 258
 259 //=======================================================================
 260 //function : Interpolate
 261 //purpose  :
 262 //=======================================================================
 263
 264 void GeomAPI_PointsToBSplineSurface::Interpolate(const TColgp_Array2OfPnt& Points,
 265                                                  const Standard_Boolean thePeriodic)
 266 {
 267   Interpolate(Points, Approx_ChordLength, thePeriodic);
 268 }
 269
 270 //=======================================================================
 271 //function : Interpolate
 272 //purpose  :
 273 //=======================================================================
 274
 275 void GeomAPI_PointsToBSplineSurface::Interpolate(const TColgp_Array2OfPnt& Points,
 276                                           const Approx_ParametrizationType ParType,
 277                                           const Standard_Boolean thePeriodic)
 278 {
 279   Standard_Integer DegMin, DegMax;
 280   DegMin = DegMax = 3;
 281   GeomAbs_Shape CC  = GeomAbs_C2;
 282   Standard_Real Tol3d = -1.0;
 283   Init(Points, ParType, DegMin, DegMax, CC, Tol3d, thePeriodic);
 284 }
 285
 286
 287 //=======================================================================
 288 //function : Init
 289 //purpose  :
 290 //=======================================================================
 291
 292 void GeomAPI_PointsToBSplineSurface::Init(const TColgp_Array2OfPnt& Points,
 293                                           const Standard_Integer DegMin,
 294                                           const Standard_Integer DegMax,
 295                                           const GeomAbs_Shape Continuity,
 296                                           const Standard_Real Tol3D)
 297 {
 298   Init(Points, Approx_ChordLength, DegMin, DegMax, Continuity, Tol3D);
 299 }
 300 //=======================================================================
 301 //function : Init
 302 //purpose  :
 303 //=======================================================================
 304
 305 void GeomAPI_PointsToBSplineSurface::Init(const TColgp_Array2OfPnt& Points,
 306             const Approx_ParametrizationType ParType,
 307             const Standard_Integer DegMin,
 308             const Standard_Integer DegMax,
 309             const GeomAbs_Shape Continuity,
 310             const Standard_Real Tol3D,
 311             const Standard_Boolean thePeriodic)
 312 {
 313   Standard_Integer Imin = Points.LowerRow();
 314   Standard_Integer Imax = Points.UpperRow();
 315   Standard_Integer Jmin = Points.LowerCol();
 316   Standard_Integer Jmax = Points.UpperCol();
 317
 318   Standard_Real Tol2D = Tol3D;
 319
 320   // first approximate the U isos:
 321   Standard_Integer add = 1;
 322   if (thePeriodic)
 323   {
 324     add = 2;
 325   }
 326   AppDef_MultiLine Line(Jmax-Jmin+1);
 327   Standard_Integer i, j;
 328
 329   for (j = Jmin; j <= Jmax; j++) {
 330     AppDef_MultiPointConstraint MP(Imax-Imin+add, 0);
 331     for (i = Imin; i <= Imax; i++) {
 332       MP.SetPoint(i, Points(i,j));
 333     }
 334     if (thePeriodic)
 335     {
 336       MP.SetPoint(Imax+1, Points(1, j));
 337     }
 338     Line.SetValue(j, MP);
 339   }
 340
 341
 342
 343
 344   Standard_Integer nbit = 2;
 345   Standard_Boolean UseSquares = Standard_False;
 346   if(Tol3D <= 1.e-3) UseSquares = Standard_True;
 347
 348   AppDef_BSplineCompute TheComputer
 349     (DegMin,DegMax,Tol3D,Tol2D,nbit,Standard_True,ParType,UseSquares);
 350
 351   switch( Continuity) {
 352   case GeomAbs_C0:
 353     TheComputer.SetContinuity(0); break;
 354
 355   case GeomAbs_G1:
 356   case GeomAbs_C1:
 357     TheComputer.SetContinuity(1); break;
 358
 359   case GeomAbs_G2:
 360   case GeomAbs_C2:
 361     TheComputer.SetContinuity(2); break;
 362
 363   default:
 364     TheComputer.SetContinuity(3);
 365   }
 366
 367   if (Tol3D <= 0.0) {
 368     TheComputer.Interpol(Line);
 369   }
 370   else {
 371     TheComputer.Perform(Line);
 372   }
 373
 374   const AppParCurves_MultiBSpCurve& TheCurve = TheComputer.Value();
 375
 376   Standard_Integer VDegree = TheCurve.Degree();
 377   TColgp_Array1OfPnt Poles(1,TheCurve.NbPoles());
 378   const TColStd_Array1OfReal& VKnots = TheCurve.Knots();
 379   const TColStd_Array1OfInteger& VMults = TheCurve.Multiplicities();
 380
 381
 382   Standard_Integer nbisosu = Imax-Imin+add;
 383   AppDef_MultiLine Line2(nbisosu);
 384
 385   for (i = 1; i <= nbisosu; i++) {
 386     TheCurve.Curve(i, Poles);
 387     AppDef_MultiPointConstraint MP(Poles.Upper(),0);
 388     for (j = 1; j <= Poles.Upper(); j++) {
 389       MP.SetPoint(j, Poles(j));
 390     }
 391     Line2.SetValue(i, MP);
 392   }
 393
 394
 395   // approximate the resulting poles:
 396
 397   AppDef_BSplineCompute TheComputer2
 398     (DegMin,DegMax,Tol3D,Tol2D,nbit,Standard_True,ParType,UseSquares);
 399   if (Tol3D <= 0.0) {
 400     if (thePeriodic)
 401     {
 402       TheComputer2.SetPeriodic(thePeriodic);
 403     }
 404     TheComputer2.Interpol(Line2);
 405   }
 406   else {
 407     if (thePeriodic && Line2.NbMultiPoints() > 2)
 408     {
 409       TheComputer2.SetPeriodic(thePeriodic);
 410       //
 411       TColStd_Array1OfReal aPars(1, Line2.NbMultiPoints());
 412       BuildParameters(Line2, ParType, aPars);
 413       math_Vector aTang(1, 3 * Poles.Upper());
 414       BuildPeriodicTangent(Line2, aPars, aTang);
 415       Standard_Integer ind = 1;
 416       TheCurve.Curve(ind, Poles);
 417       AppDef_MultiPointConstraint MP1(Poles.Upper(), 0);
 418       for (j = 1; j <= Poles.Upper(); j++) {
 419         MP1.SetPoint(j, Poles(j));
 420         Standard_Integer k = 3 * (j - 1);
 421         gp_Vec aT(aTang(k + 1), aTang(k + 2), aTang(k + 3));
 422         MP1.SetTang(j, aT);
 423       }
 424       Line2.SetValue(ind, MP1);
 425       //
 426       ind = Line2.NbMultiPoints();
 427       TheCurve.Curve(ind, Poles);
 428       AppDef_MultiPointConstraint MP2(Poles.Upper(), 0);
 429       for (j = 1; j <= Poles.Upper(); j++) {
 430         MP2.SetPoint(j, Poles(j));
 431         Standard_Integer k = 3 * (j - 1);
 432         gp_Vec aT(aTang(k + 1), aTang(k + 2), aTang(k + 3));
 433         MP2.SetTang(j, aT);
 434       }
 435       Line2.SetValue(ind, MP2);
 436     }
 437     TheComputer2.Perform(Line2);
 438   }
 439
 440   const AppParCurves_MultiBSpCurve& TheCurve2 = TheComputer2.Value();
 441
 442   Standard_Integer UDegree = TheCurve2.Degree();
 443   TColgp_Array1OfPnt Poles2(1,TheCurve2.NbPoles());
 444   const TColStd_Array1OfReal& UKnots = TheCurve2.Knots();
 445   const TColStd_Array1OfInteger& UMults = TheCurve2.Multiplicities();
 446
 447
 448   // computing the surface
 449   TColgp_Array2OfPnt ThePoles(1, Poles2.Upper(), 1, Poles.Upper());
 450   for ( j = 1; j <= Poles.Upper(); j++) {
 451     TheCurve2.Curve(j, Poles2);
 452     for (i = 1; i<= Poles2.Upper(); i++) {
 453       ThePoles(i, j) = Poles2(i);
 454     }
 455   }
 456
 457
 458   mySurface = new Geom_BSplineSurface(ThePoles, UKnots, VKnots, UMults, VMults,
 459                                       UDegree, VDegree);
 460   if (thePeriodic && Line2.NbMultiPoints() > 2)
 461   {
 462     mySurface->SetUPeriodic();
 463   }
 464
 465   myIsDone = Standard_True;
 466 }
 467
 468
 469 //=======================================================================
 470 //function : Init
 471 //purpose  :
 472 //=======================================================================
 473
 474 void GeomAPI_PointsToBSplineSurface::Init(const TColgp_Array2OfPnt& Points,
 475                                           const Standard_Real Weight1,
 476                                           const Standard_Real Weight2,
 477                                           const Standard_Real Weight3,
 478                                           const Standard_Integer DegMax,
 479                                           const GeomAbs_Shape Continuity,
 480                                           const Standard_Real Tol3D)
 481 {
 482   Standard_Integer Imin = Points.LowerRow();
 483   Standard_Integer Imax = Points.UpperRow();
 484   Standard_Integer Jmin = Points.LowerCol();
 485   Standard_Integer Jmax = Points.UpperCol();
 486
 487   Standard_Integer nbit = 2;
 488   if(Tol3D <= 1.e-3) nbit = 0;
 489
 490   // first approximate the U isos:
 491   Standard_Integer NbPointJ = Jmax-Jmin+1;
 492   Standard_Integer NbPointI = Imax-Imin+1;
 493   Standard_Integer i, j;
 494
 495   AppDef_MultiLine Line(NbPointJ);
 496
 497   for (j = Jmin; j <= Jmax; j++) {
 498     AppDef_MultiPointConstraint MP(NbPointI, 0);
 499     for (i = Imin; i <= Imax; i++) {
 500       MP.SetPoint(i, Points(i,j));
 501     }
 502     Line.SetValue(j, MP);
 503   }
 504
 505
 506   Handle(AppParCurves_HArray1OfConstraintCouple) TABofCC =
 507     new AppParCurves_HArray1OfConstraintCouple(1, NbPointJ);
 508   AppParCurves_Constraint  Constraint=AppParCurves_NoConstraint;
 509
 510   for(i = 1; i <= NbPointJ; ++i) {
 511     AppParCurves_ConstraintCouple ACC(i,Constraint);
 512     TABofCC->SetValue(i,ACC);
 513   }
 514
 515
 516   AppDef_Variational Variation(Line, 1, NbPointJ, TABofCC);
 517
 518 //===================================
 519   Standard_Integer theMaxSegments = 1000;
 520   Standard_Boolean theWithMinMax = Standard_False;
 521 //===================================
 522
 523   Variation.SetMaxDegree(DegMax);
 524   Variation.SetContinuity(Continuity);
 525   Variation.SetMaxSegment(theMaxSegments);
 526
 527   Variation.SetTolerance(Tol3D);
 528   Variation.SetWithMinMax(theWithMinMax);
 529   Variation.SetNbIterations(nbit);
 530
 531   Variation.SetCriteriumWeight(Weight1, Weight2, Weight3);
 532
 533   if(!Variation.IsCreated()) {
 534     return;
 535   }
 536
 537   if(Variation.IsOverConstrained()) {
 538     return;
 539   }
 540
 541   try {
 542     Variation.Approximate();
 543   }
 544   catch (Standard_Failure const&) {
 545     return;
 546   }
 547
 548   if(!Variation.IsDone()) {
 549     return;
 550   }
 551
 552   const AppParCurves_MultiBSpCurve& TheCurve = Variation.Value();
 553
 554   Standard_Integer VDegree = TheCurve.Degree();
 555   TColgp_Array1OfPnt Poles(1,TheCurve.NbPoles());
 556   const TColStd_Array1OfReal& VKnots = TheCurve.Knots();
 557   const TColStd_Array1OfInteger& VMults = TheCurve.Multiplicities();
 558
 559
 560   AppDef_MultiLine Line2(NbPointI);
 561
 562   for (i = 1; i <= NbPointI; i++) {
 563     TheCurve.Curve(i, Poles);
 564     AppDef_MultiPointConstraint MP(Poles.Upper(),0);
 565     for (j = 1; j <= Poles.Upper(); j++) {
 566       MP.SetPoint(j, Poles(j));
 567     }
 568     Line2.SetValue(i, MP);
 569   }
 570
 571
 572   Handle(AppParCurves_HArray1OfConstraintCouple) TABofCC2 =
 573     new AppParCurves_HArray1OfConstraintCouple(1, NbPointI);
 574
 575   for(i = 1; i <= NbPointI; ++i) {
 576     AppParCurves_ConstraintCouple ACC(i,Constraint);
 577     TABofCC2->SetValue(i,ACC);
 578   }
 579
 580
 581   AppDef_Variational Variation2(Line2, 1, NbPointI, TABofCC2);
 582
 583   Variation2.SetMaxDegree(DegMax);
 584   Variation2.SetContinuity(Continuity);
 585   Variation2.SetMaxSegment(theMaxSegments);
 586
 587   Variation2.SetTolerance(Tol3D);
 588   Variation2.SetWithMinMax(theWithMinMax);
 589   Variation2.SetNbIterations(nbit);
 590
 591   Variation2.SetCriteriumWeight(Weight1, Weight2, Weight3);
 592
 593   if(!Variation2.IsCreated()) {
 594     return;
 595   }
 596
 597   if(Variation2.IsOverConstrained()) {
 598     return;
 599   }
 600
 601   try {
 602     Variation2.Approximate();
 603   }
 604   catch (Standard_Failure const&) {
 605     return;
 606   }
 607
 608   if(!Variation2.IsDone()) {
 609     return;
 610   }
 611
 612   const AppParCurves_MultiBSpCurve& TheCurve2 = Variation2.Value();
 613
 614   Standard_Integer UDegree = TheCurve2.Degree();
 615   TColgp_Array1OfPnt Poles2(1,TheCurve2.NbPoles());
 616   const TColStd_Array1OfReal& UKnots = TheCurve2.Knots();
 617   const TColStd_Array1OfInteger& UMults = TheCurve2.Multiplicities();
 618
 619
 620   // computing the surface
 621   TColgp_Array2OfPnt ThePoles(1, Poles2.Upper(), 1, Poles.Upper());
 622   for ( j = 1; j <= Poles.Upper(); j++) {
 623     TheCurve2.Curve(j, Poles2);
 624     for (i = 1; i<= Poles2.Upper(); i++) {
 625       ThePoles(i, j) = Poles2(i);
 626     }
 627   }
 628
 629
 630   mySurface = new Geom_BSplineSurface(ThePoles, UKnots, VKnots, UMults, VMults,
 631                                       UDegree, VDegree);
 632
 633   myIsDone = Standard_True;
 634 }
 635 //=======================================================================
 636 //function : Init
 637 //purpose  :
 638 //=======================================================================
 639
 640 void GeomAPI_PointsToBSplineSurface::Interpolate(const TColStd_Array2OfReal& ZPoints,
 641                                                  const Standard_Real X0,
 642                                                  const Standard_Real dX,
 643                                                  const Standard_Real Y0,
 644                                                  const Standard_Real dY)
 645 {
 646   Standard_Integer DegMin, DegMax;
 647   DegMin = DegMax = 3;
 648   Standard_Real Tol3D  = -1.0;
 649   GeomAbs_Shape CC = GeomAbs_C2;
 650   Init(ZPoints, X0, dX, Y0, dY, DegMin, DegMax, CC, Tol3D);
 651 }
 652
 653 //=======================================================================
 654 //function : Init
 655 //purpose  :
 656 //=======================================================================
 657
 658 void GeomAPI_PointsToBSplineSurface::Init(const TColStd_Array2OfReal& ZPoints,
 659                                           const Standard_Real X0,
 660                                           const Standard_Real dX,
 661                                           const Standard_Real Y0,
 662                                           const Standard_Real dY,
 663                                           const Standard_Integer DegMin,
 664                                           const Standard_Integer DegMax,
 665                                           const GeomAbs_Shape Continuity,
 666                                           const Standard_Real Tol3D)
 667 {
 668   Standard_Integer Imin = ZPoints.LowerRow();
 669   Standard_Integer Imax = ZPoints.UpperRow();
 670   Standard_Integer Jmin = ZPoints.LowerCol();
 671   Standard_Integer Jmax = ZPoints.UpperCol();
 672   Standard_Real length;
 673
 674   Standard_Real Tol2D = Tol3D;
 675
 676   // first approximate the U isos:
 677   AppDef_MultiLine Line(Jmax-Jmin+1);
 678   math_Vector Param(Jmin, Jmax);
 679   Standard_Integer i, j;
 680 //  Standard_Real X, Y;
 681   length = dY * (Jmax-Jmin);
 682
 683   for (j = Jmin; j <= Jmax; j++) {
 684     AppDef_MultiPointConstraint MP(0, Imax-Imin+1);
 685     for (i = Imin; i <= Imax; i++) {
 686       MP.SetPoint2d(i, gp_Pnt2d(0.0, ZPoints(i, j)));
 687     }
 688     Param(j) = (Standard_Real)(j-Jmin)/(Standard_Real)(Jmax-Jmin);
 689     Line.SetValue(j, MP);
 690   }
 691
 692   AppDef_BSplineCompute TheComputer
 693     (Param, DegMin,DegMax,Tol3D,Tol2D,0, Standard_True, Standard_True);
 694
 695   switch( Continuity) {
 696   case GeomAbs_C0:
 697     TheComputer.SetContinuity(0); break;
 698
 699   case GeomAbs_G1:
 700   case GeomAbs_C1:
 701     TheComputer.SetContinuity(1); break;
 702
 703   case GeomAbs_G2:
 704   case GeomAbs_C2:
 705     TheComputer.SetContinuity(2); break;
 706
 707   default:
 708     TheComputer.SetContinuity(3);
 709   }
 710
 711   if (Tol3D <= 0.0) {
 712     TheComputer.Interpol(Line);
 713   }
 714   else {
 715     TheComputer.Perform(Line);
 716   }
 717
 718
 719   const AppParCurves_MultiBSpCurve& TheCurve = TheComputer.Value();
 720
 721   Standard_Integer VDegree = TheCurve.Degree();
 722   TColgp_Array1OfPnt2d Poles(1,TheCurve.NbPoles());
 723   Standard_Integer nk = TheCurve.Knots().Length();
 724   TColStd_Array1OfReal VKnots(1,nk);
 725   TColStd_Array1OfInteger VMults(1,nk);
 726
 727
 728   // compute Y values for the poles
 729   TColStd_Array1OfReal YPoles(1,Poles.Upper());
 730
 731   // start with a line
 732   TColStd_Array1OfReal    TempPoles(1,2);
 733   TColStd_Array1OfReal    TempKnots(1,2);
 734   TColStd_Array1OfInteger TempMults(1,2);
 735   TempMults.Init(2);
 736   TempPoles(1) = Y0;
 737   TempPoles(2) = Y0 + length;
 738   TempKnots(1) = 0.;
 739   TempKnots(2) = 1.;
 740
 741   // increase the Degree
 742   TColStd_Array1OfReal    NewTempPoles(1,VDegree+1);
 743   TColStd_Array1OfReal    NewTempKnots(1,2);
 744   TColStd_Array1OfInteger NewTempMults(1,2);
 745   BSplCLib::IncreaseDegree(1,VDegree,Standard_False,1,
 746                            TempPoles,TempKnots,TempMults,
 747                            NewTempPoles,NewTempKnots,NewTempMults);
 748
 749
 750   // insert the Knots
 751   BSplCLib::InsertKnots(VDegree,Standard_False,1,
 752                         NewTempPoles,NewTempKnots,NewTempMults,
 753                         TheCurve.Knots(),&TheCurve.Multiplicities(),
 754                         YPoles,VKnots,VMults,
 755                         Epsilon(1));
 756
 757   // scale the knots
 758   for (i = 1; i <= nk; i++) {
 759     VKnots(i) = Y0 + length * VKnots(i);
 760   }
 761
 762
 763   Standard_Integer nbisosu = Imax-Imin+1;
 764   AppDef_MultiLine Line2(nbisosu);
 765   math_Vector Param2(1, nbisosu);
 766   length = dX*(Imax-Imin);
 767
 768   for (i = 1; i <= nbisosu; i++) {
 769     TheCurve.Curve(i, Poles);
 770     AppDef_MultiPointConstraint MP(0, Poles.Upper());
 771     for (j = 1; j <= Poles.Upper(); j++) {
 772       MP.SetPoint2d(j, gp_Pnt2d(0.0, Poles(j).Y()));
 773     }
 774     Param2(i) = (Standard_Real)(i-1)/(Standard_Real)(nbisosu-1);
 775     Line2.SetValue(i, MP);
 776   }
 777
 778
 779   // approximate the resulting poles:
 780
 781   AppDef_BSplineCompute TheComputer2
 782     (Param2, DegMin,DegMax,Tol3D,Tol2D,0, Standard_True, Standard_True);
 783   if (Tol3D <= 0.0) {
 784     TheComputer2.Interpol(Line2);
 785   }
 786   else {
 787     TheComputer2.Perform(Line2);
 788   }
 789
 790   const AppParCurves_MultiBSpCurve& TheCurve2 = TheComputer2.Value();
 791
 792   Standard_Integer UDegree = TheCurve2.Degree();
 793   TColgp_Array1OfPnt2d Poles2(1,TheCurve2.NbPoles());
 794   Standard_Integer nk2 = TheCurve2.Knots().Length();
 795   TColStd_Array1OfReal UKnots(1,nk2);
 796   TColStd_Array1OfInteger UMults(1,nk2);
 797
 798
 799   // compute X values for the poles
 800   TColStd_Array1OfReal XPoles(1,Poles2.Upper());
 801
 802   // start with a line
 803   TempPoles(1) = X0;
 804   TempPoles(2) = X0 + length;
 805   TempKnots(1) = 0.;
 806   TempKnots(2) = 1.;
 807   TempMults(1) = TempMults(2) = 2;
 808
 809   // increase the Degree
 810   TColStd_Array1OfReal    NewTempPoles2(1,UDegree+1);
 811   BSplCLib::IncreaseDegree(1,UDegree,Standard_False,1,
 812                            TempPoles,TempKnots,TempMults,
 813                            NewTempPoles2,NewTempKnots,NewTempMults);
 814
 815
 816   // insert the Knots
 817   BSplCLib::InsertKnots(UDegree,Standard_False,1,
 818                         NewTempPoles2,NewTempKnots,NewTempMults,
 819                         TheCurve2.Knots(),&TheCurve2.Multiplicities(),
 820                         XPoles,UKnots,UMults,
 821                         Epsilon(1));
 822
 823   // scale the knots
 824   for (i = 1; i <= nk2; i++) {
 825     UKnots(i) = X0 + length * UKnots(i);
 826   }
 827
 828   // creating the surface
 829   TColgp_Array2OfPnt ThePoles(1, Poles2.Upper(), 1, Poles.Upper());
 830
 831   for ( j = 1; j <= Poles.Upper(); j++) {
 832     TheCurve2.Curve(j, Poles2);
 833     for (i = 1; i<= Poles2.Upper(); i++) {
 834       ThePoles(i, j).SetCoord(XPoles(i), YPoles(j), Poles2(i).Y());
 835     }
 836   }
 837
 838
 839   mySurface = new Geom_BSplineSurface(ThePoles, UKnots, VKnots, UMults, VMults,
 840                                       UDegree, VDegree);
 841
 842   myIsDone = Standard_True;
 843
 844 }
 845
 846
 847 //=======================================================================
 848 //function : Handle(Geom_BSplineSurface)&
 849 //purpose  :
 850 //=======================================================================
 851
 852 const Handle(Geom_BSplineSurface)& GeomAPI_PointsToBSplineSurface::Surface()
 853 const
 854 {
 855   StdFail_NotDone_Raise_if
 856     ( !myIsDone,
 857      "GeomAPI_PointsToBSplineSurface: Surface not done");
 858
 859   return mySurface;
 860 }
 861
 862
 863
 864 //=======================================================================
 865 //function : Handle
 866 //purpose  :
 867 //=======================================================================
 868
 869 GeomAPI_PointsToBSplineSurface::operator Handle(Geom_BSplineSurface)() const
 870 {
 871   return Surface();
 872 }
 873
 874 //=======================================================================
 875 //function : IsDone
 876 //purpose  :
 877 //=======================================================================
 878
 879 Standard_Boolean GeomAPI_PointsToBSplineSurface::IsDone () const
 880 {
 881   return myIsDone;
 882 }
 883
 884