src/AppParCurves/AppParCurves_Variational.gxx

   1 // Copyright (c) 1996-1999 Matra Datavision
   2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
   3 //
   4 // The content of this file is subject to the Open CASCADE Technology Public
   5 // License Version 6.5 (the "License"). You may not use the content of this file
   6 // except in compliance with the License. Please obtain a copy of the License
   7 // at http://www.opencascade.org and read it completely before using this file.
   8 //
   9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  11 //
  12 // The Original Code and all software distributed under the License is
  13 // distributed on an "AS IS" basis, without warranty of any kind, and the
  14 // Initial Developer hereby disclaims all such warranties, including without
  15 // limitation, any warranties of merchantability, fitness for a particular
  16 // purpose or non-infringement. Please see the License for the specific terms
  17 // and conditions governing the rights and limitations under the License.
  18
  19 // Jeannine PANCIATICI le 06/06/96
  20 // Igor FEOKTISTOV 14/12/98 - correction of Approximate() and Init().
  21 // Approximation d une MultiLine de points decrite par le tool MLineTool.
  22 // avec criteres variationnels
  23
  24
  25 #define No_Standard_RangeError
  26 #define No_Standard_OutOfRange
  27 #define No_Standard_DimensionError
  28 #define No_Standard_ConstructionError
  29
  30 #include <Standard_Stream.hxx>
  31
  32 #include <AppParCurves.hxx>
  33 #include <AppParCurves_Constraint.hxx>
  34 #include <AppParCurves_HArray1OfConstraintCouple.hxx>
  35 #include <AppParCurves_Array1OfMultiPoint.hxx>
  36 #include <AppParCurves_MultiPoint.hxx>
  37 #include <AppParCurves_MultiBSpCurve.hxx>
  38 #include <Convert_CompPolynomialToPoles.hxx>
  39 #include <gp_Pnt.hxx>
  40 #include <gp_Pnt2d.hxx>
  41 #include <gp_Vec.hxx>
  42 #include <gp_Vec2d.hxx>
  43 #include <TColgp_Array1OfPnt.hxx>
  44 #include <TColgp_Array1OfPnt2d.hxx>
  45 #include <TColgp_Array1OfVec.hxx>
  46 #include <TColgp_Array1OfVec2d.hxx>
  47 #include <TColStd_Array1OfInteger.hxx>
  48 #include <TColStd_HArray1OfInteger.hxx>
  49 #include <TColStd_Array1OfReal.hxx>
  50 #include <TColStd_HArray1OfReal.hxx>
  51 #include <TColStd_HArray2OfReal.hxx>
  52 #include <StdFail_NotDone.hxx>
  53 #include <Standard_SStream.hxx>
  54 #include <Standard_NoSuchObject.hxx>
  55 #include <Precision.hxx>
  56 //#include <Smoothing.h>
  57
  58 #if defined(WNT)
  59 # include <stdio.h>
  60 # include <stdarg.h>
  61 #endif  /* WNT */
  62
  63
  64 //
  65 //=======================================================================
  66 //function : AppParCurves_Variational
  67 //purpose  : Initialization of   the   fields.
  68 //=======================================================================
  69 //
  70  AppParCurves_Variational::AppParCurves_Variational(const MultiLine& SSP,
  71                                                     const Standard_Integer FirstPoint,
  72                                                     const Standard_Integer LastPoint,
  73                                                     const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
  74                                                     const Standard_Integer MaxDegree,
  75                                                     const Standard_Integer MaxSegment,
  76                                                     const GeomAbs_Shape Continuity,
  77                                                     const Standard_Boolean WithMinMax,
  78                                                     const Standard_Boolean WithCutting,
  79                                                     const Standard_Real Tolerance,
  80                                                     const Standard_Integer NbIterations):
  81                                         mySSP(SSP),
  82                                         myFirstPoint(FirstPoint),
  83                                         myLastPoint(LastPoint),
  84                                         myConstraints(TheConstraints),
  85                                         myMaxDegree(MaxDegree),
  86                                         myMaxSegment(MaxSegment),
  87                                         myNbIterations(NbIterations),
  88                                         myTolerance(Tolerance),
  89                                         myContinuity(Continuity),
  90                                         myWithMinMax(WithMinMax),
  91                                         myWithCutting(WithCutting)
  92
  93
  94 {
  95 // Verifications:
  96     if (myMaxDegree < 1) Standard_DomainError::Raise();
  97     myMaxDegree = Min (30, myMaxDegree);
  98 //
  99     if (myMaxSegment < 1) Standard_DomainError::Raise();
 100 //
 101     if (myWithMinMax != 0 && myWithMinMax !=1 ) Standard_DomainError::Raise();
 102     if (myWithCutting != 0 && myWithCutting !=1 ) Standard_DomainError::Raise();
 103 //
 104     myIsOverConstr = Standard_False;
 105     myIsCreated    = Standard_False;
 106     myIsDone       = Standard_False;
 107     switch (myContinuity) {
 108     case GeomAbs_C0:
 109        myNivCont=0;
 110        break ;
 111     case GeomAbs_C1:
 112        myNivCont=1;
 113        break ;
 114     case GeomAbs_C2:
 115        myNivCont=2;
 116        break ;
 117     default:
 118     Standard_ConstructionError::Raise();
 119   }
 120 //
 121     myNbP2d = ToolLine::NbP2d(SSP);
 122     myNbP3d = ToolLine::NbP3d(SSP);
 123     myDimension = 2 * myNbP2d + 3* myNbP3d ;
 124 //
 125     myPercent[0]=0.4;
 126     myPercent[1]=0.2;
 127     myPercent[2]=0.4;
 128     myKnots= new TColStd_HArray1OfReal(1,2);
 129     myKnots->SetValue(1,0.);
 130     myKnots->SetValue(2,1.);
 131
 132 //  Declaration
 133 //
 134     mySmoothCriterion = new AppParCurves_MyCriterion(mySSP, myFirstPoint, myLastPoint);
 135     myParameters = new TColStd_HArray1OfReal(myFirstPoint, myLastPoint);
 136     myNbPoints=myLastPoint-myFirstPoint+1;
 137     if (myNbPoints <= 0) Standard_ConstructionError::Raise();
 138 //
 139     myTabPoints= new TColStd_HArray1OfReal(1,myDimension*myNbPoints);
 140 //
 141 //  Table of Points initialization
 142 //
 143     Standard_Integer ipoint,jp2d,jp3d,index;
 144     TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
 145     TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
 146     gp_Pnt2d P2d;
 147     gp_Pnt P3d;
 148     index=1;
 149
 150     for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++)
 151       {
 152
 153         if(myNbP2d !=0 && myNbP3d ==0 ) {
 154                           ToolLine::Value(mySSP,ipoint,TabP2d);
 155
 156                           for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++)
 157                               {  P2d = TabP2d.Value(jp2d);
 158
 159                                  myTabPoints->SetValue(index++,P2d.X());
 160                                  myTabPoints->SetValue(index++,P2d.Y());
 161                               }
 162                          }
 163         if(myNbP3d !=0 && myNbP2d == 0 ) {
 164                           ToolLine::Value(mySSP,ipoint,TabP3d);
 165
 166                           for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
 167
 168                               {  P3d=TabP3d.Value(jp3d);
 169
 170                                  myTabPoints->SetValue(index++,P3d.X());
 171                                  myTabPoints->SetValue(index++,P3d.Y());
 172                                  myTabPoints->SetValue(index++,P3d.Z());
 173
 174                               }
 175                          }
 176         if(myNbP3d !=0 && myNbP2d != 0 ) {
 177                           ToolLine::Value(mySSP,ipoint,TabP3d,TabP2d);
 178
 179                           for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
 180
 181                               {  P3d=TabP3d.Value(jp3d);
 182
 183                                  myTabPoints->SetValue(index++,P3d.X());
 184                                  myTabPoints->SetValue(index++,P3d.Y());
 185                                  myTabPoints->SetValue(index++,P3d.Z());
 186
 187                               }
 188                           for ( jp2d = 1 ; jp2d <= myNbP2d ; jp2d++)
 189
 190                               {  P2d=TabP2d.Value(jp2d);
 191
 192                                  myTabPoints->SetValue(index++,P2d.X());
 193                                  myTabPoints->SetValue(index++,P2d.Y());
 194                               }
 195                          }
 196       }
 197     Init();
 198 }
 199 //
 200 //=======================================================================
 201 //function : Init
 202 //purpose  : Initializes the tables of constraints
 203 //           and verifies if the problem is not over-constrained.
 204 //           This method is used in the Create and the method SetConstraint.
 205 //=======================================================================
 206 //
 207 void AppParCurves_Variational::Init()
 208 {
 209
 210     Standard_Integer ipoint,jp2d,jp3d,index,jndex;
 211     Standard_Integer CurMultyPoint;
 212     TColgp_Array1OfVec TabV3d(1, Max(1,myNbP3d));
 213     TColgp_Array1OfVec2d TabV2d(1, Max(1,myNbP2d));
 214     TColgp_Array1OfVec TabV3dcurv(1, Max(1,myNbP3d));
 215     TColgp_Array1OfVec2d TabV2dcurv(1, Max(1,myNbP2d));
 216
 217     gp_Vec Vt3d, Vc3d;
 218     gp_Vec2d Vt2d, Vc2d;
 219
 220     myNbConstraints=myConstraints->Length();
 221     if (myNbConstraints < 0) Standard_ConstructionError::Raise();
 222
 223     myTypConstraints = new TColStd_HArray1OfInteger(1,Max(1,2*myNbConstraints));
 224     myTabConstraints = new TColStd_HArray1OfReal(1,Max(1,2*myDimension*myNbConstraints));
 225     myTtheta = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints));
 226     myTfthet = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints));
 227
 228
 229 //
 230 // Table of types initialization
 231       Standard_Integer iconstr;
 232       index=1;
 233       jndex=1;
 234       CurMultyPoint = 1;
 235       myNbPassPoints=0;
 236       myNbTangPoints=0;
 237       myNbCurvPoints=0;
 238       AppParCurves_Constraint valcontr;
 239
 240       for ( iconstr = myConstraints->Lower() ; iconstr <= myConstraints->Upper() ; iconstr++)
 241         {
 242         ipoint=(myConstraints->Value(iconstr)).Index();
 243
 244           valcontr=(myConstraints->Value(iconstr)).Constraint();
 245           switch (valcontr) {
 246                              case AppParCurves_NoConstraint:
 247                                   CurMultyPoint -= myNbP3d * 6 + myNbP2d * 2;
 248                                   break ;
 249                              case AppParCurves_PassPoint:
 250                                   myTypConstraints->SetValue(index++,ipoint);
 251                                   myTypConstraints->SetValue(index++,0);
 252                                   myNbPassPoints++;
 253                                   if(myNbP2d !=0 ) jndex=jndex+4*myNbP2d;
 254                                   if(myNbP3d !=0 ) jndex=jndex+6*myNbP3d;
 255                                   break ;
 256                              case AppParCurves_TangencyPoint:
 257                                   myTypConstraints->SetValue(index++,ipoint);
 258                                   myTypConstraints->SetValue(index++,1);
 259                                   myNbTangPoints++;
 260                                   if(myNbP2d !=0 && myNbP3d == 0 )
 261                                     {
 262                                      if (ToolLine::Tangency(mySSP,ipoint,TabV2d) == Standard_False)
 263                                           Standard_ConstructionError::Raise();
 264                                      for (jp2d=1;jp2d<=myNbP2d;jp2d++)
 265                                          {
 266                                           Vt2d=TabV2d.Value(jp2d);
 267                                           Vt2d.Normalize();
 268                                           myTabConstraints->SetValue(jndex++,Vt2d.X());
 269                                           myTabConstraints->SetValue(jndex++,Vt2d.Y());
 270                                           jndex=jndex+2;
 271                                           InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4);
 272                                          }
 273                                      }
 274                                    if(myNbP3d !=0 && myNbP2d == 0)
 275                                      {
 276                                       if (ToolLine::Tangency(mySSP,ipoint,TabV3d) == Standard_False)
 277                                           Standard_ConstructionError::Raise();
 278                                       for (jp3d=1;jp3d<=myNbP3d;jp3d++)
 279                                           {
 280                                            Vt3d=TabV3d.Value(jp3d);
 281                                            Vt3d.Normalize();
 282                                            myTabConstraints->SetValue(jndex++,Vt3d.X());
 283
 284                                            myTabConstraints->SetValue(jndex++,Vt3d.Y());
 285
 286                                            myTabConstraints->SetValue(jndex++,Vt3d.Z());
 287                                            jndex=jndex+3;
 288                                            InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
 289                                           }
 290                                       }
 291                                    if(myNbP3d !=0 && myNbP2d != 0)
 292                                      {
 293                                       if (ToolLine::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False)
 294                                           Standard_ConstructionError::Raise();
 295                                       for (jp3d=1;jp3d<=myNbP3d;jp3d++)
 296                                           {
 297                                            Vt3d=TabV3d.Value(jp3d);
 298                                            Vt3d.Normalize();
 299                                            myTabConstraints->SetValue(jndex++,Vt3d.X());
 300                                            myTabConstraints->SetValue(jndex++,Vt3d.Y());
 301                                            myTabConstraints->SetValue(jndex++,Vt3d.Z());
 302                                            jndex=jndex+3;
 303                                            InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
 304                                           }
 305
 306                                      for (jp2d=1;jp2d<=myNbP2d;jp2d++)
 307                                          {
 308                                           Vt2d=TabV2d.Value(jp2d);
 309                                           Vt2d.Normalize();
 310                                           myTabConstraints->SetValue(jndex++,Vt2d.X());
 311                                           myTabConstraints->SetValue(jndex++,Vt2d.Y());
 312                                           jndex=jndex+2;
 313                                           InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4);
 314                                          }
 315                                      }
 316
 317
 318                                        break ;
 319
 320                              case AppParCurves_CurvaturePoint:
 321                                    myTypConstraints->SetValue(index++,ipoint);
 322                                    myTypConstraints->SetValue(index++,2);
 323                                    myNbCurvPoints++;
 324                                    if(myNbP2d !=0 && myNbP3d == 0)
 325                                        {
 326                                        if (ToolLine::Tangency(mySSP,ipoint,TabV2d) == Standard_False )
 327                                            Standard_ConstructionError::Raise();
 328                                        if (ToolLine::Curvature(mySSP,ipoint,TabV2dcurv) == Standard_False)
 329                                            Standard_ConstructionError::Raise();
 330                                        for (jp2d=1;jp2d<=myNbP2d;jp2d++)
 331                                            {
 332                                             Vt2d=TabV2d.Value(jp2d);
 333                                             Vt2d.Normalize();
 334                                             Vc2d=TabV2dcurv.Value(jp2d);
 335                                             if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular())
 336                                                  Standard_ConstructionError::Raise();
 337                                             myTabConstraints->SetValue(jndex++,Vt2d.X());
 338                                             myTabConstraints->SetValue(jndex++,Vt2d.Y());
 339                                             myTabConstraints->SetValue(jndex++,Vc2d.X());
 340                                             myTabConstraints->SetValue(jndex++,Vc2d.Y());
 341                                             InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4);
 342                                            }
 343                                        }
 344
 345                                     if(myNbP3d !=0 && myNbP2d == 0 )
 346                                        {
 347                                        if (ToolLine::Tangency(mySSP,ipoint,TabV3d) == Standard_False )
 348                                            Standard_ConstructionError::Raise();
 349                                        if (ToolLine::Curvature(mySSP,ipoint,TabV3dcurv) == Standard_False)
 350                                            Standard_ConstructionError::Raise();
 351                                        for (jp3d=1;jp3d<=myNbP3d;jp3d++)
 352                                            {
 353                                             Vt3d=TabV3d.Value(jp3d);
 354                                             Vt3d.Normalize();
 355                                             Vc3d=TabV3dcurv.Value(jp3d);
 356                                             if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False)
 357                                                  Standard_ConstructionError::Raise();
 358                                             myTabConstraints->SetValue(jndex++,Vt3d.X());
 359                                             myTabConstraints->SetValue(jndex++,Vt3d.Y());
 360                                             myTabConstraints->SetValue(jndex++,Vt3d.Z());
 361                                             myTabConstraints->SetValue(jndex++,Vc3d.X());
 362                                             myTabConstraints->SetValue(jndex++,Vc3d.Y());
 363                                             myTabConstraints->SetValue(jndex++,Vc3d.Z());
 364                                             InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
 365                                            }
 366                                        }
 367                                     if(myNbP3d !=0 && myNbP2d != 0 )
 368                                        {
 369                                        if (ToolLine::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False )
 370                                            Standard_ConstructionError::Raise();
 371                                        if (ToolLine::Curvature(mySSP,ipoint,TabV3dcurv,TabV2dcurv) == Standard_False)
 372                                            Standard_ConstructionError::Raise();
 373                                        for (jp3d=1;jp3d<=myNbP3d;jp3d++)
 374                                            {
 375                                             Vt3d=TabV3d.Value(jp3d);
 376                                             Vt3d.Normalize();
 377                                             Vc3d=TabV3dcurv.Value(jp3d);
 378                                             if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False)
 379                                                  Standard_ConstructionError::Raise();
 380                                             myTabConstraints->SetValue(jndex++,Vt3d.X());
 381                                             myTabConstraints->SetValue(jndex++,Vt3d.Y());
 382                                             myTabConstraints->SetValue(jndex++,Vt3d.Z());
 383                                             myTabConstraints->SetValue(jndex++,Vc3d.X());
 384                                             myTabConstraints->SetValue(jndex++,Vc3d.Y());
 385                                             myTabConstraints->SetValue(jndex++,Vc3d.Z());
 386                                             InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
 387                                            }
 388                                        for (jp2d=1;jp2d<=myNbP2d;jp2d++)
 389                                            {
 390                                             Vt2d=TabV2d.Value(jp2d);
 391                                             Vt2d.Normalize();
 392                                             Vc2d=TabV2dcurv.Value(jp2d);
 393                                             if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular())
 394                                                  Standard_ConstructionError::Raise();
 395                                             myTabConstraints->SetValue(jndex++,Vt2d.X());
 396                                             myTabConstraints->SetValue(jndex++,Vt2d.Y());
 397                                             myTabConstraints->SetValue(jndex++,Vc2d.X());
 398                                             myTabConstraints->SetValue(jndex++,Vc2d.Y());
 399                                             InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4);
 400                                            }
 401
 402                                       }
 403                                       break ;
 404                                  default:
 405                                       Standard_ConstructionError::Raise();
 406                                 }
 407           CurMultyPoint += myNbP3d * 6 + myNbP2d * 2;
 408            }
 409 // OverConstraint Detection
 410         Standard_Integer MaxSeg;
 411         if(myWithCutting == Standard_True) MaxSeg = myMaxSegment ;
 412         else MaxSeg = 1;
 413         if (((myMaxDegree-myNivCont)*MaxSeg-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
 414              {
 415                myIsOverConstr =Standard_True;
 416                myIsCreated = Standard_False;
 417              }
 418         else {
 419           InitSmoothCriterion();
 420           myIsCreated = Standard_True;
 421         }
 422
 423
 424 }
 425 //
 426 //=======================================================================
 427 //function : Approximate
 428 //purpose  : Makes the approximation with the current fields.
 429 //=======================================================================
 430 //
 431 void AppParCurves_Variational::Approximate()
 432
 433 {
 434   if (myIsCreated == Standard_False )  StdFail_NotDone:: Raise();
 435
 436
 437   Standard_Real WQuadratic, WQuality;
 438
 439   TColStd_Array1OfReal Ecarts(myFirstPoint, myLastPoint);
 440
 441   mySmoothCriterion->GetWeight(WQuadratic, WQuality);
 442
 443   Handle(FEmTool_Curve) TheCurve;
 444
 445   mySmoothCriterion->GetCurve(TheCurve);
 446
 447 //---------------------------------------------------------------------
 448
 449   TheMotor(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts);
 450
 451
 452   if(myWithMinMax && myTolerance < myMaxError)
 453     Adjusting(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts);
 454
 455 //---------------------------------------------------------------------
 456
 457   Standard_Integer jp2d,jp3d,index,ipole,
 458                    NbElem = TheCurve->NbElements();
 459
 460   TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
 461   TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
 462   Standard_Real debfin[2] = {-1., 1};
 463
 464   gp_Pnt2d P2d;
 465   gp_Pnt P3d;
 466
 467   index=0;
 468
 469   {
 470     Handle(TColStd_HArray2OfReal) PolynomialIntervalsPtr =
 471       new TColStd_HArray2OfReal(1,NbElem,1,2) ;
 472
 473     Handle(TColStd_HArray1OfInteger) NbCoeffPtr =
 474       new TColStd_HArray1OfInteger(1, myMaxSegment);
 475
 476     Standard_Integer size = myMaxSegment * (myMaxDegree + 1) * myDimension ;
 477     Handle(TColStd_HArray1OfReal) CoeffPtr = new TColStd_HArray1OfReal(1,size);
 478
 479     CoeffPtr->Init(0.);
 480
 481     Handle(TColStd_HArray1OfReal) IntervallesPtr =
 482       new TColStd_HArray1OfReal(1, NbElem + 1);
 483
 484     IntervallesPtr->ChangeArray1() = TheCurve->Knots();
 485
 486     TheCurve->GetPolynom(CoeffPtr->ChangeArray1());
 487
 488     Standard_Integer ii;
 489
 490     for(ii = 1; ii <= NbElem; ii++)
 491       NbCoeffPtr->SetValue(ii, TheCurve->Degree(ii)+1);
 492
 493
 494     for (ii = PolynomialIntervalsPtr->LowerRow() ;
 495          ii <= PolynomialIntervalsPtr->UpperRow() ;ii++)
 496       {
 497         PolynomialIntervalsPtr->SetValue(ii,1,debfin[0]) ;
 498         PolynomialIntervalsPtr->SetValue(ii,2,debfin[1]) ;
 499       }
 500 /*
 501         printf("\n =========== Parameters for converting\n");
 502         printf("nb_courbes : %d \n", NbElem);
 503         printf("tab_option[4] : %d \n", myNivCont);
 504         printf("myDimension : %d \n", myDimension);
 505         printf("myMaxDegree : %d \n", myMaxDegree);
 506         printf("\n NbcoeffPtr :\n");
 507         for(ii = 1; ii <= NbElem; ii++)  printf("NbCoeffPtr(%d) = %d \n", ii, NbCoeffPtr->Value(ii));
 508         size = NbElem*(myMaxDegree + 1) * myDimension;
 509         printf("\n CoeffPtr :\n");
 510         for(ii = 1; ii <= size; ii++)  printf("CoeffPtr(%d) = %.8e \n", ii, CoeffPtr->Value(ii));
 511         printf("\n PolinimialIntervalsPtr :\n");
 512         for (ii = PolynomialIntervalsPtr->LowerRow() ;
 513              ii <= PolynomialIntervalsPtr->UpperRow() ;ii++)
 514             {
 515              printf(" %d : [%f, %f] \n", ii, PolynomialIntervalsPtr->Value(ii,1),
 516                     PolynomialIntervalsPtr->Value(ii,2)) ;
 517             }
 518         printf("\n IntervallesPtr :\n");
 519         for (ii = IntervallesPtr->Lower();
 520              ii <= IntervallesPtr->Upper() - 1; ii++)
 521             {
 522              printf(" %d : [%f, %f] \n", ii, IntervallesPtr->Value(ii),
 523                     IntervallesPtr->Value(ii+1)) ;
 524             }
 525 */
 526     Convert_CompPolynomialToPoles AConverter(NbElem,
 527                                              myNivCont,
 528                                              myDimension,
 529                                              myMaxDegree,
 530                                              NbCoeffPtr,
 531                                              CoeffPtr,
 532                                              PolynomialIntervalsPtr,
 533                                              IntervallesPtr) ;
 534     if (AConverter.IsDone())
 535       {
 536         Handle(TColStd_HArray2OfReal) PolesPtr ;
 537         Handle(TColStd_HArray1OfInteger) Mults;
 538         Standard_Integer NbPoles=AConverter.NbPoles();
 539 //      Standard_Integer Deg=AConverter.Degree();
 540         AppParCurves_Array1OfMultiPoint TabMU(1,NbPoles);
 541         AConverter.Poles(PolesPtr) ;
 542         AConverter.Knots(myKnots) ;
 543         AConverter.Multiplicities(Mults) ;
 544
 545         for (ipole=PolesPtr->LowerRow();ipole<=PolesPtr->UpperRow();ipole++)
 546           {
 547             index=PolesPtr->LowerCol();
 548 /*          if(myNbP2d !=0 )
 549               {
 550                 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
 551                   {
 552                     P2d.SetX(PolesPtr->Value(ipole,index++));
 553                     P2d.SetY(PolesPtr->Value(ipole,index++));
 554                     TabP2d.SetValue(jp2d,P2d);
 555                   }
 556               }*/
 557             if(myNbP3d !=0 )
 558               {
 559                 for (jp3d=1;jp3d<=myNbP3d;jp3d++)
 560                   {
 561                     //                       cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
 562                     P3d.SetX(PolesPtr->Value(ipole,index++));
 563                     //                       cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
 564                     P3d.SetY(PolesPtr->Value(ipole,index++));
 565                     //                       cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
 566                     P3d.SetZ(PolesPtr->Value(ipole,index++));
 567                     TabP3d.SetValue(jp3d,P3d);
 568                   }
 569               }
 570             if(myNbP2d !=0 )
 571               {
 572                 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
 573                   {
 574                     P2d.SetX(PolesPtr->Value(ipole,index++));
 575                     P2d.SetY(PolesPtr->Value(ipole,index++));
 576                     TabP2d.SetValue(jp2d,P2d);
 577                   }
 578               }
 579             if(myNbP2d !=0 && myNbP3d !=0)
 580               {
 581                 AppParCurves_MultiPoint aMultiPoint(TabP3d,TabP2d);
 582                 TabMU.SetValue(ipole,aMultiPoint);
 583               }
 584             else if (myNbP2d !=0)
 585               {
 586                 AppParCurves_MultiPoint aMultiPoint(TabP2d);
 587                 TabMU.SetValue(ipole,aMultiPoint);
 588               }
 589             else
 590               {
 591
 592                 AppParCurves_MultiPoint aMultiPoint(TabP3d);
 593                 TabMU.SetValue(ipole,aMultiPoint);
 594               }
 595
 596           }
 597         AppParCurves_MultiBSpCurve aCurve(TabMU,myKnots->Array1(),Mults->Array1());
 598         myMBSpCurve=aCurve;
 599         myIsDone = Standard_True;
 600       }
 601   }
 602
 603 }
 604 //
 605 //=======================================================================
 606 //function : IsCreated
 607 //purpose  : returns True if the creation is done
 608 //=======================================================================
 609 //
 610 Standard_Boolean AppParCurves_Variational::IsCreated() const
 611 {
 612   return myIsCreated;
 613 }
 614 //
 615 //=======================================================================
 616 //function : IsDone
 617 //purpose  : returns True if the  approximation is ok
 618 //=======================================================================
 619 //
 620 Standard_Boolean AppParCurves_Variational::IsDone() const
 621 {
 622   return myIsDone;
 623 }
 624 //
 625 //=======================================================================
 626 //function : IsOverConstrained
 627 //purpose  : returns True if the problem is overconstrained
 628 //           in this case, approximation cannot be done.
 629 //=======================================================================
 630 //
 631 Standard_Boolean AppParCurves_Variational::IsOverConstrained() const
 632 {
 633   return myIsOverConstr;
 634 }
 635 //
 636 //=======================================================================
 637 //function : Value
 638 //purpose  : returns all the BSpline curves approximating the
 639 //           MultiLine SSP after minimization of the parameter.
 640
 641 //=======================================================================
 642 //
 643 AppParCurves_MultiBSpCurve AppParCurves_Variational::Value() const
 644 {
 645   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
 646   return myMBSpCurve;
 647
 648 }
 649 //
 650 //=======================================================================
 651 //function : MaxError
 652 //purpose  : returns the maximum of the distances between
 653 //           the points of the multiline and the approximation
 654 //           curves.
 655 //=======================================================================
 656 //
 657 Standard_Real AppParCurves_Variational::MaxError() const
 658 {
 659   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
 660   return myMaxError;
 661 }
 662 //
 663 //=======================================================================
 664 //function : MaxErrorIndex
 665 //purpose  : returns the index of the MultiPoint of ErrorMax
 666 //=======================================================================
 667 //
 668 Standard_Integer AppParCurves_Variational::MaxErrorIndex() const
 669 {
 670   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
 671   return myMaxErrorIndex;
 672 }
 673 //
 674 //=======================================================================
 675 //function : QuadraticError
 676 //purpose  :  returns the quadratic average of the distances between
 677 //            the points of the multiline and the approximation
 678 //            curves.
 679 //=======================================================================
 680 //
 681 Standard_Real AppParCurves_Variational::QuadraticError() const
 682 {
 683   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
 684   return myCriterium[0];
 685 }
 686 //
 687 //=======================================================================
 688 //function : Distance
 689 //purpose  : returns the distances between the points of the
 690 //           multiline and the approximation curves.
 691 //=======================================================================
 692 //
 693 void AppParCurves_Variational::Distance(math_Matrix& mat)
 694
 695 {
 696       if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
 697       Standard_Integer ipoint,jp2d,jp3d,index;
 698       TColgp_Array1OfPnt TabP3d(1,Max(1,myNbP3d));
 699       TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
 700       Standard_Integer j0 = mat.LowerCol() - myFirstPoint;
 701
 702       gp_Pnt2d P2d;
 703       gp_Pnt P3d;
 704
 705
 706       gp_Pnt Pt3d;
 707       gp_Pnt2d Pt2d;
 708
 709       for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++)
 710         {
 711           index=1;
 712           if(myNbP3d !=0 ) {
 713                             ToolLine::Value(mySSP,ipoint,TabP3d);
 714
 715                             for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
 716
 717                                       {  P3d=TabP3d.Value(jp3d);
 718                                          myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt3d);
 719                                          mat(index++, j0 + ipoint)=P3d.Distance(Pt3d);
 720
 721                                       }
 722                             }
 723           if(myNbP2d !=0 ) {
 724                             if(myNbP3d == 0 ) ToolLine::Value(mySSP,ipoint,TabP2d);
 725                             else ToolLine::Value(mySSP,ipoint,TabP3d,TabP2d);
 726                             for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++)
 727
 728                                       {  P2d = TabP2d.Value(jp2d);
 729                                          myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt2d);
 730                                          mat(index++, j0 + ipoint)=P2d.Distance(Pt2d);
 731                                       }
 732                            }
 733       }
 734
 735 }
 736 //
 737 //=======================================================================
 738 //function : AverageError
 739 //purpose  : returns the average error between
 740 //           the MultiLine and the approximation.
 741 //=======================================================================
 742 //
 743 Standard_Real AppParCurves_Variational::AverageError() const
 744 {
 745   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
 746   return myAverageError;
 747 }
 748 //
 749 //=======================================================================
 750 //function : Parameters
 751 //purpose  : returns the parameters uses to the approximations
 752 //=======================================================================
 753 //
 754 const Handle(TColStd_HArray1OfReal)& AppParCurves_Variational::Parameters() const
 755 {
 756   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
 757   return myParameters;
 758 }
 759 //
 760 //=======================================================================
 761 //function : Knots
 762 //purpose  : returns the knots uses to the approximations
 763 //=======================================================================
 764 //
 765 const Handle(TColStd_HArray1OfReal)& AppParCurves_Variational::Knots() const
 766 {
 767   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
 768   return myKnots;
 769 }
 770 //
 771 //=======================================================================
 772 //function : Criterium
 773 //purpose  : returns the values of the quality criterium.
 774 //=======================================================================
 775 //
 776 void AppParCurves_Variational::Criterium(Standard_Real& VFirstOrder, Standard_Real& VSecondOrder, Standard_Real& VThirdOrder) const
 777 {
 778   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
 779   VFirstOrder=myCriterium[1] ;
 780   VSecondOrder=myCriterium[2];
 781   VThirdOrder=myCriterium[3];
 782 }
 783 //
 784 //=======================================================================
 785 //function : CriteriumWeight
 786 //purpose  : returns the Weights (as percent) associed  to the criterium used in
 787 //           the  optimization.
 788 //=======================================================================
 789 //
 790 void AppParCurves_Variational::CriteriumWeight(Standard_Real& Percent1, Standard_Real& Percent2, Standard_Real& Percent3) const
 791 {
 792   Percent1 = myPercent[0];
 793   Percent2 = myPercent[1];
 794   Percent3 = myPercent[2] ;
 795 }
 796 //
 797 //=======================================================================
 798 //function : MaxDegree
 799 //purpose  : returns the Maximum Degree used in the approximation
 800 //=======================================================================
 801 //
 802 Standard_Integer AppParCurves_Variational::MaxDegree() const
 803 {
 804    return myMaxDegree;
 805 }
 806 //
 807 //=======================================================================
 808 //function : MaxSegment
 809 //purpose  : returns the Maximum of segment used in the approximation
 810 //=======================================================================
 811 //
 812 Standard_Integer AppParCurves_Variational::MaxSegment() const
 813 {
 814     return myMaxSegment;
 815  }
 816 //
 817 //=======================================================================
 818 //function : Continuity
 819 //purpose  : returns the Continuity used in the approximation
 820 //=======================================================================
 821 //
 822 GeomAbs_Shape AppParCurves_Variational::Continuity() const
 823 {
 824     return myContinuity;
 825 }
 826 //
 827 //=======================================================================
 828 //function : WithMinMax
 829 //purpose  : returns if the  approximation  search to  minimize the
 830 //           maximum Error or not.
 831 //=======================================================================
 832 //
 833 Standard_Boolean AppParCurves_Variational::WithMinMax() const
 834 {
 835     return myWithMinMax;
 836 }
 837 //
 838 //=======================================================================
 839 //function : WithCutting
 840 //purpose  :  returns if the  approximation can insert new Knots or not.
 841 //=======================================================================
 842 //
 843 Standard_Boolean AppParCurves_Variational::WithCutting() const
 844 {
 845    return myWithCutting;
 846 }
 847 //
 848 //=======================================================================
 849 //function : Tolerance
 850 //purpose  : returns the tolerance used in the approximation.
 851 //=======================================================================
 852 //
 853 Standard_Real AppParCurves_Variational::Tolerance() const
 854 {
 855    return myTolerance;
 856 }
 857 //
 858 //=======================================================================
 859 //function : NbIterations
 860 //purpose  : returns the number of iterations used in the approximation.
 861 //=======================================================================
 862 //
 863 Standard_Integer AppParCurves_Variational::NbIterations() const
 864 {
 865    return myNbIterations;
 866 }
 867 //
 868 //=======================================================================
 869 //function : Dump
 870 //purpose  : Prints on the stream o information on the current state
 871 //           of the object.
 872 //=======================================================================
 873 //
 874 void AppParCurves_Variational::Dump(Standard_OStream& o) const
 875 {
 876   o << " \nVariational Smoothing " << endl;
 877   o << " Number of multipoints                   "  << myNbPoints << endl;
 878   o << " Number of 2d par multipoint "  << myNbP2d << endl;
 879   o << " Nombre of 3d par multipoint "  << myNbP3d << endl;
 880   o << " Number of PassagePoint      "  << myNbPassPoints << endl;
 881   o << " Number of TangencyPoints    "  << myNbTangPoints << endl;
 882   o << " Number of CurvaturePoints   "  << myNbCurvPoints << endl;
 883   o << " \nTolerance " << o.setf(ios::scientific) << setprecision(3) << setw(9) << myTolerance;
 884   if ( WithMinMax()) { o << "  as Max Error." << endl;}
 885   else { o << "  as size Error." << endl;}
 886   o << "CriteriumWeights : " << myPercent[0] << " , "
 887     << myPercent[1] << " , " << myPercent[2] << endl;
 888
 889   if (myIsDone ) {
 890     o << " MaxError             " << setprecision(3) << setw(9) << myMaxError << endl;
 891     o << " Index of  MaxError   " << myMaxErrorIndex << endl;
 892     o << " Average Error        " << setprecision(3) << setw(9) << myAverageError << endl;
 893     o << " Quadratic Error      " << setprecision(3) << setw(9) << myCriterium[0] << endl;
 894     o << " Tension Criterium    " << setprecision(3) << setw(9) << myCriterium[1] << endl;
 895     o << " Flexion  Criterium   " << setprecision(3) << setw(9) << myCriterium[2] << endl;
 896     o << " Jerk  Criterium      " << setprecision(3) << setw(9) << myCriterium[3] << endl;
 897     o << " NbSegments           "  << myKnots->Length()-1 << endl;
 898   }
 899   else
 900    { if (myIsOverConstr) o << "The probleme is overconstraint " << endl;
 901      else o << " Erreur dans l''approximation" << endl;
 902    }
 903 }
 904 //
 905 //=======================================================================
 906 //function : SetConstraints
 907 //purpose  : Define the constraints to approximate
 908 //           If this value is incompatible with the others fields
 909 //           this method modify nothing and returns false
 910 //=======================================================================
 911 //
 912 Standard_Boolean AppParCurves_Variational::SetConstraints(const Handle(AppParCurves_HArray1OfConstraintCouple)& aConstraint)
 913
 914 {
 915    myConstraints=aConstraint;
 916    Init();
 917    if (myIsOverConstr ) return Standard_False;
 918    else return Standard_True;
 919 }
 920 //
 921 //=======================================================================
 922 //function : SetParameters
 923 //purpose  : Defines the parameters used by the approximations.
 924 //=======================================================================
 925 //
 926 void AppParCurves_Variational::SetParameters(const Handle(TColStd_HArray1OfReal)& param)
 927 {
 928   myParameters->ChangeArray1() = param->Array1();
 929 }
 930 //
 931 //=======================================================================
 932 //function : SetKnots
 933 //purpose  : Defines the knots used by the approximations
 934 //    --          If this value is incompatible with the others fields
 935 //    --          this method modify nothing and returns false
 936 //=======================================================================
 937 //
 938 Standard_Boolean AppParCurves_Variational::SetKnots(const Handle(TColStd_HArray1OfReal)& knots)
 939 {
 940   myKnots->ChangeArray1() = knots->Array1();
 941   return Standard_True;
 942 }
 943 //
 944 //=======================================================================
 945 //function : SetMaxDegree
 946 //purpose  : Define the Maximum Degree used in the approximation
 947 //           If this value is incompatible with the others fields
 948 //           this method modify nothing and returns false
 949 //=======================================================================
 950 //
 951 Standard_Boolean AppParCurves_Variational::SetMaxDegree(const Standard_Integer Degree)
 952 {
 953   if (((Degree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
 954        return Standard_False;
 955   else
 956      {
 957           myMaxDegree=Degree;
 958
 959           InitSmoothCriterion();
 960
 961           return Standard_True;
 962         }
 963
 964 }
 965 //
 966 //=======================================================================
 967 //function : SetMaxSegment
 968 //purpose  : Define the maximum number of segments used in the approximation
 969 //           If this value is incompatible with the others fields
 970 //           this method modify nothing and returns false
 971 //=======================================================================
 972 //
 973 Standard_Boolean AppParCurves_Variational::SetMaxSegment(const Standard_Integer NbSegment)
 974 {
 975   if ( myWithCutting == Standard_True &&
 976       ((myMaxDegree-myNivCont)*NbSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
 977        return Standard_False;
 978   else
 979     {
 980        myMaxSegment=NbSegment;
 981        return Standard_True;
 982     }
 983 }
 984 //
 985 //=======================================================================
 986 //function : SetContinuity
 987 //purpose  : Define the Continuity used in the approximation
 988 //           If this value is incompatible with the others fields
 989 //           this method modify nothing and returns false
 990 //=======================================================================
 991 //
 992 Standard_Boolean AppParCurves_Variational::SetContinuity(const GeomAbs_Shape C)
 993 {
 994     Standard_Integer NivCont=0;
 995     switch (C) {
 996     case GeomAbs_C0:
 997        NivCont=0;
 998        break ;
 999     case GeomAbs_C1:
1000        NivCont=1;
1001        break ;
1002     case GeomAbs_C2:
1003        NivCont=2;
1004        break ;
1005     default:
1006     Standard_ConstructionError::Raise();
1007   }
1008   if (((myMaxDegree-NivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1009        return Standard_False;
1010   else
1011     {
1012          myContinuity= C;
1013          myNivCont=NivCont;
1014
1015          InitSmoothCriterion();
1016        return Standard_True;
1017     }
1018 }
1019 //
1020 //=======================================================================
1021 //function : SetWithMinMax
1022 //purpose  : Define if the  approximation  search to  minimize the
1023 //           maximum Error or not.
1024 //=======================================================================
1025 //
1026 void AppParCurves_Variational::SetWithMinMax(const Standard_Boolean MinMax)
1027 {
1028   myWithMinMax=MinMax;
1029
1030   InitSmoothCriterion();
1031 }
1032 //
1033 //=======================================================================
1034 //function : SetWithCutting
1035 //purpose  : Define if the  approximation can insert new Knots or not.
1036 //           If this value is incompatible with the others fields
1037 //           this method modify nothing and returns false
1038 //=======================================================================
1039 //
1040 Standard_Boolean AppParCurves_Variational::SetWithCutting(const Standard_Boolean Cutting)
1041 {
1042   if (Cutting == Standard_False)
1043       {
1044        if (((myMaxDegree-myNivCont)*myKnots->Length()-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1045        return Standard_False;
1046
1047        else
1048           {
1049         myWithCutting=Cutting;
1050         InitSmoothCriterion();
1051         return Standard_True;
1052       }
1053       }
1054   else
1055     {
1056       if (((myMaxDegree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1057        return Standard_False;
1058
1059        else
1060       {
1061         myWithCutting=Cutting;
1062         InitSmoothCriterion();
1063         return Standard_True;
1064       }
1065     }
1066 }
1067 //
1068 //=======================================================================
1069 //function : SetCriteriumWeight
1070 //purpose  : define the Weights (as percent) associed to the criterium used in
1071 //           the  optimization.
1072 //=======================================================================
1073 //
1074 void AppParCurves_Variational::SetCriteriumWeight(const Standard_Real Percent1, const Standard_Real Percent2, const Standard_Real Percent3)
1075 {
1076   if (Percent1 < 0 || Percent2 < 0 || Percent3 < 0 ) Standard_DomainError::Raise();
1077   Standard_Real Total = Percent1 + Percent2 + Percent3;
1078   myPercent[0] = Percent1/Total;
1079   myPercent[1] = Percent2/Total;
1080   myPercent[2] = Percent3/Total;
1081
1082   InitSmoothCriterion();
1083 }
1084 //
1085 //=======================================================================
1086 //function : SetCriteriumWeight
1087 //purpose  :  define the  Weight   (as  percent)  associed  to   the
1088 //            criterium   Order used in   the optimization  : Others
1089 //            weights are updated.
1090 //=======================================================================
1091 //
1092 void AppParCurves_Variational::SetCriteriumWeight(const Standard_Integer Order, const Standard_Real Percent)
1093 {
1094   if ( Percent < 0 ) Standard_DomainError::Raise();
1095   if ( Order < 1 || Order > 3 ) Standard_ConstructionError::Raise();
1096   myPercent[Order-1] = Percent;
1097   Standard_Real Total = myPercent[0] + myPercent[1] + myPercent[2];
1098   myPercent[0] = myPercent[0] / Total;
1099   myPercent[1] = myPercent[1] / Total;
1100   myPercent[2] = myPercent[2] / Total;
1101
1102   InitSmoothCriterion();
1103
1104 }
1105 //
1106 //=======================================================================
1107 //function : SetTolerance
1108 //purpose  : define the tolerance used in the approximation.
1109 //=======================================================================
1110 //
1111 void AppParCurves_Variational::SetTolerance(const Standard_Real Tol)
1112 {
1113   myTolerance=Tol;
1114   InitSmoothCriterion();
1115 }
1116 //
1117 //=======================================================================
1118 //function : SetNbIterations
1119 //purpose  : define the number of iterations used in the approximation.
1120 //=======================================================================
1121 //
1122 void AppParCurves_Variational::SetNbIterations(const Standard_Integer Iter)
1123 {
1124   myNbIterations=Iter;
1125 }
1126
1127
1128 // Private methods
1129 #include <AppParCurves_Variational_1.gxx> // TheMotor
1130 #include <AppParCurves_Variational_2.gxx> // Optimization
1131 #include <AppParCurves_Variational_3.gxx> // Project
1132 #include <AppParCurves_Variational_4.gxx> // ACR
1133 #include <AppParCurves_Variational_5.gxx> // SplitCurve
1134 #include <AppParCurves_Variational_6.gxx> // InitSmoothCriterion and other Init... methods
1135 #include <AppParCurves_Variational_7.gxx> // Adjusting
1136 #include <AppParCurves_Variational_8.gxx> // AssemblingConstraints
1137 #include <AppParCurves_Variational_9.gxx> // InitTtheta and InitTfthet methods