src/GeomPlate/GeomPlate_BuildAveragePlane.cxx

   1 // Created on: 1996-03-18
   2 // Created by: Stagiaire Frederic CALOONE
   3 // Copyright (c) 1996-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 // Modified:    Wed Mar  5 09:45:42 1997
  18 //    by:       Joelle CHAUVET
  19 //              ajout de la methode DefPlan et des options POption et NOption
  20 //              modif des methodes Create et BasePlan
  21 // Modified:    Thu Mar 20 09:15:36 1997
  22 //    by:       Joelle CHAUVET
  23 //              correction sur le tri des valeurs propres quand valeurs egales
  24
  25
  26 #include <GeomPlate_BuildAveragePlane.ixx>
  27 #include <TColgp_Array1OfVec.hxx>
  28 #include <math_Matrix.hxx>
  29 #include <gp_Pnt.hxx>
  30 #include <gp_Pln.hxx>
  31 #include <TColgp_HArray1OfPnt.hxx>
  32 #include <math_Jacobi.hxx>
  33 #include <gp_Vec.hxx>
  34 #include <gp_Ax1.hxx>
  35 #include <gp_Ax3.hxx>
  36 #include <gp_Dir.hxx>
  37 #include <ElSLib.hxx>
  38 #include <Geom_Line.hxx>
  39 #include <TColgp_Array1OfPnt2d.hxx>
  40 #include <TColStd_Array1OfReal.hxx>
  41 #include <gp_Lin2d.hxx>
  42 #include <ElCLib.hxx>
  43 #include <GeomLib.hxx>
  44
  45 #include <GeomPlate_Aij.hxx>
  46
  47
  48 //=======================================================================
  49 //function : GeomPlate_BuildAveragePlane
  50 //purpose  :
  51 //=======================================================================
  52
  53
  54 GeomPlate_BuildAveragePlane::
  55 GeomPlate_BuildAveragePlane(const Handle(TColgp_HArray1OfPnt)& Pts,
  56                             const Standard_Integer NbBoundPoints,
  57                             const Standard_Real Tol,
  58                             const Standard_Integer POption,
  59                             const Standard_Integer NOption) :
  60 myPts(Pts),
  61 myTol(Tol),
  62 myNbBoundPoints(NbBoundPoints)
  63 {
  64   gp_Vec OZ = DefPlan(NOption);
  65
  66   if (OZ.SquareMagnitude()>0) {
  67
  68     if (POption==1) {
  69       myPlane = new Geom_Plane(myG,OZ);
  70       myOX = myPlane->Pln().XAxis().Direction();
  71       myOY = myPlane->Pln().YAxis().Direction();
  72     }
  73     else {
  74       BasePlan(OZ);
  75       gp_Dir NDir(myOX^myOY);
  76       gp_Dir UDir(myOX);
  77       gp_Ax3 triedre(myG,NDir,UDir);
  78       myPlane = new Geom_Plane(triedre);
  79     }
  80     Standard_Integer i,nb=myPts->Length();
  81     gp_Pln P=myPlane->Pln();
  82     ElSLib::Parameters(P,myG,myUmax,myVmax);
  83     myUmin=myUmax;
  84     myVmin=myVmax;
  85     Standard_Real U=0,V=0;
  86     for (i=1; i<=nb; i++) {
  87       ElSLib::Parameters(P,myPts->Value(i),U,V);
  88       if ( myUmax < U ) myUmax=U;
  89       if ( myUmin > U ) myUmin=U;
  90       if ( myVmax < V ) myVmax=V;
  91       if ( myVmin > V ) myVmin=V;
  92     }
  93   }
  94
  95   if (IsLine()) {
  96     myLine = new Geom_Line(myG,myOX);
  97   }
  98 }
  99
 100 GeomPlate_BuildAveragePlane::GeomPlate_BuildAveragePlane( const TColgp_SequenceOfVec& Normals,
 101                                                           const Handle( TColgp_HArray1OfPnt )& Pts ) :
 102 myPts(Pts)
 103 {
 104   Standard_Integer i, j, k, n, m;
 105
 106   gp_Vec BestVec;
 107   Standard_Integer NN = Normals.Length();
 108
 109   if (NN == 1)
 110     BestVec = Normals(1);
 111   else if (NN == 2)
 112     {
 113       BestVec = Normals(1) + Normals(2);
 114       BestVec.Normalize();
 115     }
 116   else //the common case
 117     {
 118       Standard_Real MaxAngle = 0.;
 119       for (i = 1; i <= NN-1; i++)
 120         for (j = i+1; j <= NN; j++)
 121           {
 122             Standard_Real Angle = Normals(i).Angle( Normals(j) );
 123             if (Angle > MaxAngle)
 124               MaxAngle = Angle;
 125           }
 126       MaxAngle *= 1.2;
 127       MaxAngle /= 2.;
 128       Standard_Integer Nint = 50;
 129
 130       TColgp_Array1OfVec OptVec( 1, NN*(NN-1)/2 );
 131       TColStd_Array1OfReal OptScal( 1, NN*(NN-1)/2 );
 132       gp_Vec Vec, Vec1;
 133       gp_Dir Cross1, Cross2;
 134
 135       k = 1;
 136       for (i = 1; i <= NN-1; i++)
 137         for (j = i+1; j <= NN; j++)
 138           {
 139             Standard_Real Step = MaxAngle/Nint;
 140             Vec = Normals(i) + Normals(j);
 141             Vec.Normalize();
 142
 143             Cross1 = Normals(i) ^ Normals(j);
 144             Cross2 = Vec ^ Cross1;
 145             gp_Ax1 Axe( gp_Pnt(0,0,0), Cross2 );
 146
 147             Vec1 = Vec.Rotated( Axe, -MaxAngle );
 148             //Vec2 = Vec.Rotated( Axe, MaxAngle );
 149
 150             OptScal(k) = RealFirst();
 151             for (n = 0; n <= 2*Nint; n++)
 152               {
 153                 Vec1.Rotate( Axe, Step );
 154                 Standard_Real minScal = RealLast();
 155                 for (m = 1; m <= NN; m++)
 156                   {
 157                     Standard_Real Scal = Vec1 * Normals(m);
 158                     if (Scal < minScal)
 159                       minScal = Scal;
 160                   }
 161                 if (minScal > OptScal(k))
 162                   {
 163                     OptScal(k) = minScal;
 164                     OptVec(k) = Vec1;
 165                   }
 166               }
 167             k++;
 168           } // for i, for j
 169       //Find maximum among all maximums
 170
 171       Standard_Real BestScal = RealFirst();
 172       Standard_Integer Index=0;
 173       for (k = 1; k <= OptScal.Length(); k++)
 174         if (OptScal(k) > BestScal)
 175           {
 176             BestScal = OptScal(k);
 177             Index = k;
 178           }
 179       BestVec = OptVec(Index);
 180     }
 181
 182   //Making the plane myPlane
 183   gp_Ax2 Axe;
 184   Standard_Boolean IsSingular;
 185   TColgp_Array1OfPnt PtsArray( 1, myPts->Length() );
 186   for (i = 1; i <= myPts->Length(); i++)
 187     PtsArray(i) = myPts->Value(i);
 188   GeomLib::AxeOfInertia( PtsArray, Axe, IsSingular );
 189   gp_Dir BestDir( BestVec );
 190   gp_Dir XDir = BestDir ^ Axe.XDirection();
 191   XDir ^= BestDir;
 192
 193   gp_Ax3 Axe3( Axe.Location(), BestDir, XDir );
 194   myPlane = new Geom_Plane( Axe3 );
 195
 196   //Initializing myUmin, myVmin, myUmax, myVmax
 197   gp_Pln Pln = myPlane->Pln();
 198   ElSLib::Parameters( Pln, Axe.Location(), myUmax, myVmax );
 199   myUmin = myUmax;
 200   myVmin = myVmax;
 201   Standard_Real U,V;
 202   for (i = 1; i <= myPts->Length(); i++)
 203     {
 204       gp_Vec aVec( Pln.Location(), myPts->Value(i) );
 205       gp_Vec NormVec = Pln.Axis().Direction();
 206       NormVec = (aVec * NormVec) * NormVec;
 207
 208       ElSLib::Parameters( Pln, myPts->Value(i).Translated( -NormVec ), U, V ); //????? Real projecting?
 209       if (U > myUmax)
 210         myUmax = U;
 211       if (U < myUmin)
 212         myUmin = U;
 213       if (V > myVmax)
 214         myVmax = V;
 215       if (V < myVmin)
 216         myVmin = V;
 217     }
 218   //Initializing myOX, myOY
 219   myOX = myPlane->Pln().XAxis().Direction();
 220   myOY = myPlane->Pln().YAxis().Direction();
 221 }
 222
 223
 224 //=======================================================================
 225 //function : Plane
 226 //purpose  :
 227 //=======================================================================
 228
 229 Handle(Geom_Plane) GeomPlate_BuildAveragePlane::Plane() const
 230 {
 231   Standard_NoSuchObject_Raise_if ( IsLine() , "Cannot use the function 'GeomPlate_BuildAveragePlane::Plane()', the Object is a 'Geom_Line'");
 232   return myPlane;
 233 }
 234
 235
 236 //=======================================================================
 237 //function : MinMaxBox
 238 //purpose  :
 239 //=======================================================================
 240
 241 void GeomPlate_BuildAveragePlane::MinMaxBox(Standard_Real& Umin, Standard_Real& Umax, Standard_Real& Vmin, Standard_Real& Vmax) const
 242 {
 243   Umax=myUmax;
 244   Umin=myUmin;
 245   Vmax=myVmax;
 246   Vmin=myVmin;
 247 }
 248
 249
 250
 251
 252
 253 //=======================================================================
 254 //function : DefPlan
 255 //purpose  :
 256 //=======================================================================
 257
 258 gp_Vec GeomPlate_BuildAveragePlane::DefPlan(const Standard_Integer NOption)
 259 {
 260
 261   gp_Pnt GB;
 262   gp_Vec A,B,C,D;
 263   gp_Vec OZ;
 264   Standard_Integer i,nb=myPts->Length();
 265   GB.SetCoord(0.,0.,0.);
 266   for (i=1; i<=nb; i++) {
 267       GB.SetCoord(1,(GB.Coord(1)+myPts->Value(i).Coord(1)));
 268       GB.SetCoord(2,(GB.Coord(2)+myPts->Value(i).Coord(2)));
 269       GB.SetCoord(3,(GB.Coord(3)+myPts->Value(i).Coord(3)));
 270   }
 271   myG.SetCoord(1,(GB.Coord(1)/nb));
 272   myG.SetCoord(2,(GB.Coord(2)/nb));
 273   myG.SetCoord(3,(GB.Coord(3)/nb));
 274
 275   if (NOption==1) {
 276     gp_Ax2 Axe;
 277     Standard_Boolean IsSingular;
 278     GeomLib::AxeOfInertia( myPts->Array1(), Axe, IsSingular, myTol );
 279
 280     myOX = Axe.XDirection();
 281     myOY = Axe.YDirection();
 282
 283     OZ = Axe.Direction();
 284
 285     if (myNbBoundPoints != 0 && myPts->Length() != myNbBoundPoints)
 286       {
 287         A.SetCoord(0.,0.,0.);
 288         for (i = 3; i <= myNbBoundPoints; i++)
 289           {
 290             B.SetCoord(myPts->Value(i-1).Coord(1)-myPts->Value(1).Coord(1),
 291                        myPts->Value(i-1).Coord(2)-myPts->Value(1).Coord(2),
 292                        myPts->Value(i-1).Coord(3)-myPts->Value(1).Coord(3));
 293             C.SetCoord(myPts->Value(i).Coord(1)-myPts->Value(1).Coord(1),
 294                        myPts->Value(i).Coord(2)-myPts->Value(1).Coord(2),
 295                        myPts->Value(i).Coord(3)-myPts->Value(1).Coord(3));
 296             D=B^C;
 297             A.SetCoord(1,A.Coord(1)+D.Coord(1));
 298             A.SetCoord(2,A.Coord(2)+D.Coord(2));
 299             A.SetCoord(3,A.Coord(3)+D.Coord(3));
 300           }
 301         gp_Vec OZ1 = A;
 302         Standard_Real theAngle = OZ.Angle( OZ1 );
 303         if (theAngle > M_PI/2)
 304           theAngle = M_PI - theAngle;
 305         if (theAngle > M_PI/3)
 306           OZ = OZ1;
 307       }
 308   }
 309
 310   else if (NOption==2) {
 311     A.SetCoord(0.,0.,0.);
 312     for (i = 3; i <= myNbBoundPoints; i++) {
 313       B.SetCoord(myPts->Value(i-1).Coord(1)-myPts->Value(1).Coord(1),
 314                  myPts->Value(i-1).Coord(2)-myPts->Value(1).Coord(2),
 315                  myPts->Value(i-1).Coord(3)-myPts->Value(1).Coord(3));
 316       C.SetCoord(myPts->Value(i).Coord(1)-myPts->Value(1).Coord(1),
 317                  myPts->Value(i).Coord(2)-myPts->Value(1).Coord(2),
 318                  myPts->Value(i).Coord(3)-myPts->Value(1).Coord(3));
 319       D=B^C;
 320       A.SetCoord(1,A.Coord(1)+D.Coord(1));
 321       A.SetCoord(2,A.Coord(2)+D.Coord(2));
 322       A.SetCoord(3,A.Coord(3)+D.Coord(3));
 323     }
 324     OZ = A;
 325   }
 326   return OZ;
 327 }
 328
 329 //=======================================================================
 330 //function : BasePlan
 331 //purpose  :
 332 //=======================================================================
 333
 334 void GeomPlate_BuildAveragePlane::BasePlan(const gp_Vec& OZ)
 335 {
 336     math_Matrix M (1, 3, 1, 3);
 337     M.Init(0.);
 338     gp_Vec Proj;
 339     Standard_Integer i,nb=myPts->Length();
 340     Standard_Real scal;
 341
 342     for (i=1; i<=nb; i++) {
 343       Proj.SetCoord(1,myPts->Value(i).Coord(1) - myG.Coord(1));
 344       Proj.SetCoord(2,myPts->Value(i).Coord(2) - myG.Coord(2));
 345       Proj.SetCoord(3,myPts->Value(i).Coord(3) - myG.Coord(3));
 346       scal = Proj.Coord(1)*OZ.Coord(1)
 347         + Proj.Coord(2)*OZ.Coord(2)
 348         + Proj.Coord(3)*OZ.Coord(3);
 349       scal /= OZ.Coord(1)*OZ.Coord(1)
 350         + OZ.Coord(2)*OZ.Coord(2)
 351         + OZ.Coord(3)*OZ.Coord(3);
 352       Proj.SetCoord(1,Proj.Coord(1) - scal*OZ.Coord(1));
 353       Proj.SetCoord(2,Proj.Coord(2) - scal*OZ.Coord(2));
 354       Proj.SetCoord(3,Proj.Coord(3) - scal*OZ.Coord(3));
 355       M(1,1) += Proj.Coord(1)*Proj.Coord(1);
 356       M(2,2) += Proj.Coord(2)*Proj.Coord(2);
 357       M(3,3) += Proj.Coord(3)*Proj.Coord(3);
 358       M(1,2) += Proj.Coord(1)*Proj.Coord(2);
 359       M(1,3) += Proj.Coord(1)*Proj.Coord(3);
 360       M(2,3) += Proj.Coord(2)*Proj.Coord(3);
 361     }
 362     M(2,1) = M(1,2) ;
 363     M(3,1) = M(1,3) ;
 364     M(3,2) = M(2,3) ;
 365     math_Jacobi J(M);
 366     Standard_Real n1,n2,n3;
 367     math_Vector V1(1,3),V2(1,3),V3(1,3);
 368     n1=J.Value(1);
 369     n2=J.Value(2);
 370     n3=J.Value(3);
 371
 372     Standard_Real r1 = Min(Min(n1,n2),n3), r2;
 373     Standard_Integer m1, m2, m3;
 374     if (r1==n1) {
 375       m1 = 1;
 376       r2 = Min(n2,n3);
 377       if (r2==n2) {
 378         m2 = 2;
 379         m3 = 3;
 380       }
 381       else {
 382         m2 = 3;
 383         m3 = 2;
 384       }
 385     }
 386     else {
 387       if (r1==n2) {
 388         m1 = 2 ;
 389         r2 = Min(n1,n3);
 390         if (r2==n1) {
 391           m2 = 1;
 392           m3 = 3;
 393         }
 394         else {
 395           m2 = 3;
 396           m3 = 1;
 397         }
 398       }
 399       else {
 400         m1 = 3 ;
 401         r2 = Min(n1,n2);
 402         if (r2==n1) {
 403           m2 = 1;
 404           m3 = 2;
 405         }
 406         else {
 407           m2 = 2;
 408           m3 = 1;
 409         }
 410       }
 411     }
 412     J.Vector(m1,V1);
 413     J.Vector(m2,V2);
 414     J.Vector(m3,V3);
 415
 416     if (((Abs(n1)<=myTol)&&(Abs(n2)<=myTol))
 417           || ((Abs(n2)<=myTol)&&(Abs(n3)<=myTol))
 418           || ((Abs(n1)<=myTol)&&(Abs(n3)<=myTol))) {
 419       myOX.SetCoord(V3(1),V3(2),V3(3));
 420       myOY.SetCoord(0,0,0);
 421     }
 422     else {
 423       myOX.SetCoord(V3(1),V3(2),V3(3));
 424       myOY.SetCoord(V2(1),V2(2),V2(3));
 425     }
 426 }
 427
 428
 429
 430 //=======================================================================
 431 //function : Line
 432 //purpose  :
 433 //=======================================================================
 434
 435 Handle(Geom_Line) GeomPlate_BuildAveragePlane::Line() const
 436 {
 437   Standard_NoSuchObject_Raise_if ( IsPlane() , "Cannot use the function 'GeomPlate_BuildAveragePlane::Line()', the Object is a 'Geom_Plane'");
 438   return myLine;
 439 }
 440
 441 //=======================================================================
 442 //function : IsPlane
 443 //purpose  :
 444 //=======================================================================
 445
 446 Standard_Boolean GeomPlate_BuildAveragePlane::IsPlane() const
 447 {
 448   gp_Vec OZ=myOX^myOY;
 449   if (OZ.SquareMagnitude()==0)
 450     return Standard_False;
 451   else
 452     return Standard_True;
 453
 454 }
 455
 456 //=======================================================================
 457 //function : IsLine
 458 //purpose  :
 459 //=======================================================================
 460
 461 Standard_Boolean GeomPlate_BuildAveragePlane::IsLine() const
 462 {
 463   gp_Vec OZ=myOX^myOY;
 464   if (OZ.SquareMagnitude()==0)
 465     return Standard_True;
 466   else
 467     return Standard_False;
 468 }
 469
 470
 471 Standard_Boolean GeomPlate_BuildAveragePlane::HalfSpace( const TColgp_SequenceOfVec& NewNormals,
 472                                                          TColgp_SequenceOfVec& Normals,
 473                                                          GeomPlate_SequenceOfAij& Bset,
 474                                                          const Standard_Real LinTol,
 475                                                          const Standard_Real AngTol )
 476 {
 477   Standard_Real SquareTol = LinTol * LinTol;
 478
 479   TColgp_SequenceOfVec SaveNormals;
 480   GeomPlate_SequenceOfAij SaveBset;
 481   // 1
 482   SaveNormals = Normals;
 483   SaveBset = Bset;
 484
 485   gp_Vec Cross, NullVec( 0, 0, 0 );
 486   GeomPlate_SequenceOfAij B1set, B2set;
 487   Standard_Integer i, j, k;
 488
 489   i = 1;
 490   if (Normals.IsEmpty())
 491     {
 492       if (NewNormals.Length() == 1)
 493         {
 494           Normals.Append( NewNormals.Last() );
 495           return Standard_True;
 496         }
 497       // 2
 498       Cross = NewNormals(1) ^ NewNormals(2);
 499       if (Cross.SquareMagnitude() <= SquareTol)
 500         return Standard_False;
 501
 502       Cross.Normalize();
 503       Bset.Append( GeomPlate_Aij( 1, 2, Cross ) );
 504       Bset.Append( GeomPlate_Aij( 2, 1, -Cross) );
 505       Normals.Append( NewNormals(1) );
 506       Normals.Append( NewNormals(2) );
 507       i = 3;
 508     }
 509
 510   for (; i <= NewNormals.Length(); i++)
 511     {
 512       // 3
 513       Standard_Real Scal;
 514       for (j = 1; j <= Bset.Length(); j++)
 515         if ((Scal = Bset(j).Vec * NewNormals(i)) >= -LinTol)
 516           B2set.Append( Bset(j) );
 517
 518       Standard_Integer ii = Normals.Length()+1;
 519       for (j = 1; j <= ii-1; j++)
 520         {
 521           if (Normals(j).SquareMagnitude() == 0.)
 522             continue;
 523           // 4
 524           Cross = NewNormals(i) ^ Normals(j);
 525           if (Cross.SquareMagnitude() <= SquareTol)
 526             {
 527               Normals = SaveNormals;
 528               Bset = SaveBset;
 529               return Standard_False;
 530             }
 531           Cross.Normalize();
 532           Standard_Boolean isNew = Standard_True;
 533           for (k = 1; k <= B2set.Length(); k++)
 534             if (B2set(k).Vec.IsOpposite( -Cross, AngTol )) //if (B2set(k).Vec.IsEqual( Cross, LinTol, AngTol ))
 535               {
 536                 gp_Vec Cross1, Cross2;
 537                 Standard_Integer ind1 = B2set(k).Ind1, ind2 = B2set(k).Ind2;
 538                 if (ind1 == ii || ind2 == ii)
 539                   {
 540                     isNew = Standard_False;
 541                     break;
 542                   }
 543                 Cross1 = Normals( ind1 ) ^ NewNormals(i);
 544                 Cross2 = Normals( ind2 ) ^ NewNormals(i);
 545                 if (Cross1.SquareMagnitude() <= SquareTol || Cross2.SquareMagnitude() <= SquareTol)
 546                   {
 547                     Normals = SaveNormals;
 548                     Bset = SaveBset;
 549                     return Standard_False;
 550                   }
 551                 if (Cross1.IsOpposite( Cross2, AngTol ))
 552                   {
 553                     Cross2 = Normals( ind1 ) ^ Normals( ind2 );
 554                     if (Cross1.IsOpposite( Cross2, AngTol ))
 555                       {
 556                         Normals = SaveNormals;
 557                         Bset = SaveBset;
 558                         return Standard_False;
 559                       }
 560                   }
 561                 else
 562                   {
 563                     if (NewNormals(i).Angle( Normals( ind1 ) ) > NewNormals(i).Angle( Normals( ind2 ) ))
 564                       {
 565                         B2set(k).Ind2 = ind1;
 566                         B2set(k).Ind1 = ii;
 567                       }
 568                     else
 569                       B2set(k).Ind1 = ii;
 570                   }
 571                 isNew = Standard_False;
 572                 break;
 573               }
 574           if (isNew)
 575             B1set.Append( GeomPlate_Aij( ii, j, Cross ) );
 576
 577           Cross.Reverse();
 578           isNew = Standard_True;
 579           for (k = 1; k <= B2set.Length(); k++)
 580             if (B2set(k).Vec.IsOpposite( -Cross, AngTol )) //if (B2set(k).Vec.IsEqual( Cross, LinTol, AngTol ))
 581               {
 582                 gp_Vec Cross1, Cross2;
 583                 Standard_Integer ind1 = B2set(k).Ind1, ind2 = B2set(k).Ind2;
 584                 if (ind1 == ii || ind2 == ii)
 585                   {
 586                     isNew = Standard_False;
 587                     break;
 588                   }
 589                 Cross1 = Normals( ind1 ) ^ NewNormals(i);
 590                 Cross2 = Normals( ind2 ) ^ NewNormals(i);
 591                 if (Cross1.SquareMagnitude() <= SquareTol || Cross2.SquareMagnitude() <= SquareTol)
 592                   {
 593                     Normals = SaveNormals;
 594                     Bset = SaveBset;
 595                     return Standard_False;
 596                   }
 597                 if (Cross1.IsOpposite( Cross2, AngTol ))
 598                   {
 599                     Cross2 = Normals( ind1 ) ^ Normals( ind2 );
 600                     if (Cross1.IsOpposite( Cross2, AngTol ))
 601                       {
 602                         Normals = SaveNormals;
 603                         Bset = SaveBset;
 604                         return Standard_False;
 605                       }
 606                   }
 607                 else
 608                   {
 609                     if (NewNormals(i).Angle( Normals( ind1 ) ) > NewNormals(i).Angle( Normals( ind2 ) ))
 610                       {
 611                         B2set(k).Ind2 = ind1;
 612                         B2set(k).Ind1 = ii;
 613                       }
 614                     else
 615                       B2set(k).Ind1 = ii;
 616                   }
 617                 isNew = Standard_False;
 618                 break;
 619               }
 620           if (isNew)
 621             B1set.Append( GeomPlate_Aij( ii, j, Cross ) );
 622         }
 623
 624       // 5
 625       for (j = 1; j <= B1set.Length(); j++)
 626         {
 627           Standard_Boolean isGEnull = Standard_True;
 628           for (k = 1; k <= Normals.Length(); k++)
 629             {
 630               if (Normals(k).SquareMagnitude() == 0.)
 631                 continue;
 632               if (B1set(j).Vec * Normals(k) < -LinTol)
 633                 {
 634                   isGEnull = Standard_False;
 635                   break;
 636                 }
 637             }
 638           if (isGEnull)
 639             B2set.Append( B1set(j) );
 640         }
 641
 642       // 6
 643       if (B2set.IsEmpty())
 644         {
 645           Normals = SaveNormals;
 646           Bset = SaveBset;
 647           return Standard_False;
 648         }
 649
 650       // 7
 651       Bset = B2set;
 652       B2set.Clear();
 653       B1set.Clear();
 654       Normals.Append( NewNormals(i) );
 655
 656       // 8
 657       for (j = 1; j <= Normals.Length(); j++)
 658         {
 659           if (Normals(j).SquareMagnitude() == 0.)
 660             continue;
 661           Standard_Boolean isFound = Standard_False;
 662           for (k = 1; k <= Bset.Length(); k++)
 663             if (j == Bset(k).Ind1 || j == Bset(k).Ind2)
 664               {
 665                 isFound = Standard_True;
 666                 break;
 667               }
 668           if (! isFound)
 669             Normals(j) = NullVec;
 670         }
 671     }
 672
 673   return Standard_True;
 674 }