1 // File: GeomPlate_BuildAveragePlane.cxx
2 // Created: Mon Mar 18 15:36:29 1996
3 // Author: Stagiaire Frederic CALOONE
5 // Modified: Wed Mar 5 09:45:42 1997
7 // ajout de la methode DefPlan et des options POption et NOption
8 // modif des methodes Create et BasePlan
9 // Modified: Thu Mar 20 09:15:36 1997
11 // correction sur le tri des valeurs propres quand valeurs egales
14 #include <GeomPlate_BuildAveragePlane.ixx>
15 #include <TColgp_Array1OfVec.hxx>
16 #include <math_Matrix.hxx>
19 #include <TColgp_HArray1OfPnt.hxx>
20 #include <math_Jacobi.hxx>
26 #include <Geom_Line.hxx>
27 #include <TColgp_Array1OfPnt2d.hxx>
28 #include <TColStd_Array1OfReal.hxx>
29 #include <gp_Lin2d.hxx>
31 #include <GeomLib.hxx>
33 #include <GeomPlate_Aij.hxx>
36 //=======================================================================
37 //function : GeomPlate_BuildAveragePlane
39 //=======================================================================
42 GeomPlate_BuildAveragePlane::
43 GeomPlate_BuildAveragePlane(const Handle(TColgp_HArray1OfPnt)& Pts,
44 const Standard_Integer NbBoundPoints,
45 const Standard_Real Tol,
46 const Standard_Integer POption,
47 const Standard_Integer NOption) :
50 myNbBoundPoints(NbBoundPoints)
52 gp_Vec OZ = DefPlan(NOption);
54 if (OZ.SquareMagnitude()>0) {
57 myPlane = new Geom_Plane(myG,OZ);
58 myOX = myPlane->Pln().XAxis().Direction();
59 myOY = myPlane->Pln().YAxis().Direction();
63 gp_Dir NDir(myOX^myOY);
65 gp_Ax3 triedre(myG,NDir,UDir);
66 myPlane = new Geom_Plane(triedre);
68 Standard_Integer i,nb=myPts->Length();
69 gp_Pln P=myPlane->Pln();
70 ElSLib::Parameters(P,myG,myUmax,myVmax);
73 Standard_Real U=0,V=0;
74 for (i=1; i<=nb; i++) {
75 ElSLib::Parameters(P,myPts->Value(i),U,V);
76 if ( myUmax < U ) myUmax=U;
77 if ( myUmin > U ) myUmin=U;
78 if ( myVmax < V ) myVmax=V;
79 if ( myVmin > V ) myVmin=V;
84 myLine = new Geom_Line(myG,myOX);
88 GeomPlate_BuildAveragePlane::GeomPlate_BuildAveragePlane( const TColgp_SequenceOfVec& Normals,
89 const Handle( TColgp_HArray1OfPnt )& Pts ) :
92 Standard_Integer i, j, k, n, m;
95 Standard_Integer NN = Normals.Length();
101 BestVec = Normals(1) + Normals(2);
104 else //the common case
106 Standard_Real MaxAngle = 0.;
107 for (i = 1; i <= NN-1; i++)
108 for (j = i+1; j <= NN; j++)
110 Standard_Real Angle = Normals(i).Angle( Normals(j) );
111 if (Angle > MaxAngle)
116 Standard_Integer Nint = 50;
118 TColgp_Array1OfVec OptVec( 1, NN*(NN-1)/2 );
119 TColStd_Array1OfReal OptScal( 1, NN*(NN-1)/2 );
121 gp_Dir Cross1, Cross2;
124 for (i = 1; i <= NN-1; i++)
125 for (j = i+1; j <= NN; j++)
127 Standard_Real Step = MaxAngle/Nint;
128 Vec = Normals(i) + Normals(j);
131 Cross1 = Normals(i) ^ Normals(j);
132 Cross2 = Vec ^ Cross1;
133 gp_Ax1 Axe( gp_Pnt(0,0,0), Cross2 );
135 Vec1 = Vec.Rotated( Axe, -MaxAngle );
136 //Vec2 = Vec.Rotated( Axe, MaxAngle );
138 OptScal(k) = RealFirst();
139 for (n = 0; n <= 2*Nint; n++)
141 Vec1.Rotate( Axe, Step );
142 Standard_Real minScal = RealLast();
143 for (m = 1; m <= NN; m++)
145 Standard_Real Scal = Vec1 * Normals(m);
149 if (minScal > OptScal(k))
151 OptScal(k) = minScal;
157 //Find maximum among all maximums
159 Standard_Real BestScal = RealFirst();
160 Standard_Integer Index=0;
161 for (k = 1; k <= OptScal.Length(); k++)
162 if (OptScal(k) > BestScal)
164 BestScal = OptScal(k);
167 BestVec = OptVec(Index);
170 //Making the plane myPlane
172 Standard_Boolean IsSingular;
173 TColgp_Array1OfPnt PtsArray( 1, myPts->Length() );
174 for (i = 1; i <= myPts->Length(); i++)
175 PtsArray(i) = myPts->Value(i);
176 GeomLib::AxeOfInertia( PtsArray, Axe, IsSingular );
177 gp_Dir BestDir( BestVec );
178 gp_Dir XDir = BestDir ^ Axe.XDirection();
181 gp_Ax3 Axe3( Axe.Location(), BestDir, XDir );
182 myPlane = new Geom_Plane( Axe3 );
184 //Initializing myUmin, myVmin, myUmax, myVmax
185 gp_Pln Pln = myPlane->Pln();
186 ElSLib::Parameters( Pln, Axe.Location(), myUmax, myVmax );
190 for (i = 1; i <= myPts->Length(); i++)
192 gp_Vec aVec( Pln.Location(), myPts->Value(i) );
193 gp_Vec NormVec = Pln.Axis().Direction();
194 NormVec = (aVec * NormVec) * NormVec;
196 ElSLib::Parameters( Pln, myPts->Value(i).Translated( -NormVec ), U, V ); //????? Real projecting?
206 //Initializing myOX, myOY
207 myOX = myPlane->Pln().XAxis().Direction();
208 myOY = myPlane->Pln().YAxis().Direction();
212 //=======================================================================
215 //=======================================================================
217 Handle(Geom_Plane) GeomPlate_BuildAveragePlane::Plane() const
219 Standard_NoSuchObject_Raise_if ( IsLine() , "Cannot use the function 'GeomPlate_BuildAveragePlane::Plane()', the Object is a 'Geom_Line'");
224 //=======================================================================
225 //function : MinMaxBox
227 //=======================================================================
229 void GeomPlate_BuildAveragePlane::MinMaxBox(Standard_Real& Umin, Standard_Real& Umax, Standard_Real& Vmin, Standard_Real& Vmax) const
241 //=======================================================================
244 //=======================================================================
246 gp_Vec GeomPlate_BuildAveragePlane::DefPlan(const Standard_Integer NOption)
252 Standard_Integer i,nb=myPts->Length();
253 GB.SetCoord(0.,0.,0.);
254 for (i=1; i<=nb; i++) {
255 GB.SetCoord(1,(GB.Coord(1)+myPts->Value(i).Coord(1)));
256 GB.SetCoord(2,(GB.Coord(2)+myPts->Value(i).Coord(2)));
257 GB.SetCoord(3,(GB.Coord(3)+myPts->Value(i).Coord(3)));
259 myG.SetCoord(1,(GB.Coord(1)/nb));
260 myG.SetCoord(2,(GB.Coord(2)/nb));
261 myG.SetCoord(3,(GB.Coord(3)/nb));
265 Standard_Boolean IsSingular;
266 GeomLib::AxeOfInertia( myPts->Array1(), Axe, IsSingular, myTol );
268 myOX = Axe.XDirection();
269 myOY = Axe.YDirection();
271 OZ = Axe.Direction();
273 if (myNbBoundPoints != 0 && myPts->Length() != myNbBoundPoints)
275 A.SetCoord(0.,0.,0.);
276 for (i = 3; i <= myNbBoundPoints; i++)
278 B.SetCoord(myPts->Value(i-1).Coord(1)-myPts->Value(1).Coord(1),
279 myPts->Value(i-1).Coord(2)-myPts->Value(1).Coord(2),
280 myPts->Value(i-1).Coord(3)-myPts->Value(1).Coord(3));
281 C.SetCoord(myPts->Value(i).Coord(1)-myPts->Value(1).Coord(1),
282 myPts->Value(i).Coord(2)-myPts->Value(1).Coord(2),
283 myPts->Value(i).Coord(3)-myPts->Value(1).Coord(3));
285 A.SetCoord(1,A.Coord(1)+D.Coord(1));
286 A.SetCoord(2,A.Coord(2)+D.Coord(2));
287 A.SetCoord(3,A.Coord(3)+D.Coord(3));
290 Standard_Real theAngle = OZ.Angle( OZ1 );
291 if (theAngle > M_PI/2)
292 theAngle = M_PI - theAngle;
293 if (theAngle > M_PI/3)
298 else if (NOption==2) {
299 A.SetCoord(0.,0.,0.);
300 for (i = 3; i <= myNbBoundPoints; i++) {
301 B.SetCoord(myPts->Value(i-1).Coord(1)-myPts->Value(1).Coord(1),
302 myPts->Value(i-1).Coord(2)-myPts->Value(1).Coord(2),
303 myPts->Value(i-1).Coord(3)-myPts->Value(1).Coord(3));
304 C.SetCoord(myPts->Value(i).Coord(1)-myPts->Value(1).Coord(1),
305 myPts->Value(i).Coord(2)-myPts->Value(1).Coord(2),
306 myPts->Value(i).Coord(3)-myPts->Value(1).Coord(3));
308 A.SetCoord(1,A.Coord(1)+D.Coord(1));
309 A.SetCoord(2,A.Coord(2)+D.Coord(2));
310 A.SetCoord(3,A.Coord(3)+D.Coord(3));
317 //=======================================================================
318 //function : BasePlan
320 //=======================================================================
322 void GeomPlate_BuildAveragePlane::BasePlan(const gp_Vec& OZ)
324 math_Matrix M (1, 3, 1, 3);
327 Standard_Integer i,nb=myPts->Length();
330 for (i=1; i<=nb; i++) {
331 Proj.SetCoord(1,myPts->Value(i).Coord(1) - myG.Coord(1));
332 Proj.SetCoord(2,myPts->Value(i).Coord(2) - myG.Coord(2));
333 Proj.SetCoord(3,myPts->Value(i).Coord(3) - myG.Coord(3));
334 scal = Proj.Coord(1)*OZ.Coord(1)
335 + Proj.Coord(2)*OZ.Coord(2)
336 + Proj.Coord(3)*OZ.Coord(3);
337 scal /= OZ.Coord(1)*OZ.Coord(1)
338 + OZ.Coord(2)*OZ.Coord(2)
339 + OZ.Coord(3)*OZ.Coord(3);
340 Proj.SetCoord(1,Proj.Coord(1) - scal*OZ.Coord(1));
341 Proj.SetCoord(2,Proj.Coord(2) - scal*OZ.Coord(2));
342 Proj.SetCoord(3,Proj.Coord(3) - scal*OZ.Coord(3));
343 M(1,1) += Proj.Coord(1)*Proj.Coord(1);
344 M(2,2) += Proj.Coord(2)*Proj.Coord(2);
345 M(3,3) += Proj.Coord(3)*Proj.Coord(3);
346 M(1,2) += Proj.Coord(1)*Proj.Coord(2);
347 M(1,3) += Proj.Coord(1)*Proj.Coord(3);
348 M(2,3) += Proj.Coord(2)*Proj.Coord(3);
354 Standard_Real n1,n2,n3;
355 math_Vector V1(1,3),V2(1,3),V3(1,3);
360 Standard_Real r1 = Min(Min(n1,n2),n3), r2;
361 Standard_Integer m1, m2, m3;
404 if (((Abs(n1)<=myTol)&&(Abs(n2)<=myTol))
405 || ((Abs(n2)<=myTol)&&(Abs(n3)<=myTol))
406 || ((Abs(n1)<=myTol)&&(Abs(n3)<=myTol))) {
407 myOX.SetCoord(V3(1),V3(2),V3(3));
408 myOY.SetCoord(0,0,0);
411 myOX.SetCoord(V3(1),V3(2),V3(3));
412 myOY.SetCoord(V2(1),V2(2),V2(3));
418 //=======================================================================
421 //=======================================================================
423 Handle(Geom_Line) GeomPlate_BuildAveragePlane::Line() const
425 Standard_NoSuchObject_Raise_if ( IsPlane() , "Cannot use the function 'GeomPlate_BuildAveragePlane::Line()', the Object is a 'Geom_Plane'");
429 //=======================================================================
432 //=======================================================================
434 Standard_Boolean GeomPlate_BuildAveragePlane::IsPlane() const
437 if (OZ.SquareMagnitude()==0)
438 return Standard_False;
440 return Standard_True;
444 //=======================================================================
447 //=======================================================================
449 Standard_Boolean GeomPlate_BuildAveragePlane::IsLine() const
452 if (OZ.SquareMagnitude()==0)
453 return Standard_True;
455 return Standard_False;
459 Standard_Boolean GeomPlate_BuildAveragePlane::HalfSpace( const TColgp_SequenceOfVec& NewNormals,
460 TColgp_SequenceOfVec& Normals,
461 GeomPlate_SequenceOfAij& Bset,
462 const Standard_Real LinTol,
463 const Standard_Real AngTol )
465 Standard_Real SquareTol = LinTol * LinTol;
467 TColgp_SequenceOfVec SaveNormals;
468 GeomPlate_SequenceOfAij SaveBset;
470 SaveNormals = Normals;
473 gp_Vec Cross, NullVec( 0, 0, 0 );
474 GeomPlate_SequenceOfAij B1set, B2set;
475 Standard_Integer i, j, k;
478 if (Normals.IsEmpty())
480 if (NewNormals.Length() == 1)
482 Normals.Append( NewNormals.Last() );
483 return Standard_True;
486 Cross = NewNormals(1) ^ NewNormals(2);
487 if (Cross.SquareMagnitude() <= SquareTol)
488 return Standard_False;
492 GeomPlate_Aij A1( 1, 2, Cross );
494 Bset.Append( GeomPlate_Aij( 1, 2, Cross ) );
495 Bset.Append( GeomPlate_Aij( 2, 1, -Cross) );
496 Normals.Append( NewNormals(1) );
497 Normals.Append( NewNormals(2) );
501 for (; i <= NewNormals.Length(); i++)
505 for (j = 1; j <= Bset.Length(); j++)
506 if ((Scal = Bset(j).Vec * NewNormals(i)) >= -LinTol)
507 B2set.Append( Bset(j) );
509 Standard_Integer ii = Normals.Length()+1;
510 for (j = 1; j <= ii-1; j++)
512 if (Normals(j).SquareMagnitude() == 0.)
515 Cross = NewNormals(i) ^ Normals(j);
516 if (Cross.SquareMagnitude() <= SquareTol)
518 Normals = SaveNormals;
520 return Standard_False;
523 Standard_Boolean isNew = Standard_True;
524 for (k = 1; k <= B2set.Length(); k++)
525 if (B2set(k).Vec.IsOpposite( -Cross, AngTol )) //if (B2set(k).Vec.IsEqual( Cross, LinTol, AngTol ))
527 gp_Vec Cross1, Cross2;
528 Standard_Integer ind1 = B2set(k).Ind1, ind2 = B2set(k).Ind2;
529 if (ind1 == ii || ind2 == ii)
531 isNew = Standard_False;
534 Cross1 = Normals( ind1 ) ^ NewNormals(i);
535 Cross2 = Normals( ind2 ) ^ NewNormals(i);
536 if (Cross1.SquareMagnitude() <= SquareTol || Cross2.SquareMagnitude() <= SquareTol)
538 Normals = SaveNormals;
540 return Standard_False;
542 if (Cross1.IsOpposite( Cross2, AngTol ))
544 Cross2 = Normals( ind1 ) ^ Normals( ind2 );
545 if (Cross1.IsOpposite( Cross2, AngTol ))
547 Normals = SaveNormals;
549 return Standard_False;
554 if (NewNormals(i).Angle( Normals( ind1 ) ) > NewNormals(i).Angle( Normals( ind2 ) ))
556 B2set(k).Ind2 = ind1;
562 isNew = Standard_False;
566 B1set.Append( GeomPlate_Aij( ii, j, Cross ) );
569 isNew = Standard_True;
570 for (k = 1; k <= B2set.Length(); k++)
571 if (B2set(k).Vec.IsOpposite( -Cross, AngTol )) //if (B2set(k).Vec.IsEqual( Cross, LinTol, AngTol ))
573 gp_Vec Cross1, Cross2;
574 Standard_Integer ind1 = B2set(k).Ind1, ind2 = B2set(k).Ind2;
575 if (ind1 == ii || ind2 == ii)
577 isNew = Standard_False;
580 Cross1 = Normals( ind1 ) ^ NewNormals(i);
581 Cross2 = Normals( ind2 ) ^ NewNormals(i);
582 if (Cross1.SquareMagnitude() <= SquareTol || Cross2.SquareMagnitude() <= SquareTol)
584 Normals = SaveNormals;
586 return Standard_False;
588 if (Cross1.IsOpposite( Cross2, AngTol ))
590 Cross2 = Normals( ind1 ) ^ Normals( ind2 );
591 if (Cross1.IsOpposite( Cross2, AngTol ))
593 Normals = SaveNormals;
595 return Standard_False;
600 if (NewNormals(i).Angle( Normals( ind1 ) ) > NewNormals(i).Angle( Normals( ind2 ) ))
602 B2set(k).Ind2 = ind1;
608 isNew = Standard_False;
612 B1set.Append( GeomPlate_Aij( ii, j, Cross ) );
616 for (j = 1; j <= B1set.Length(); j++)
618 Standard_Boolean isGEnull = Standard_True;
619 for (k = 1; k <= Normals.Length(); k++)
621 if (Normals(k).SquareMagnitude() == 0.)
623 if (B1set(j).Vec * Normals(k) < -LinTol)
625 isGEnull = Standard_False;
630 B2set.Append( B1set(j) );
636 Normals = SaveNormals;
638 return Standard_False;
645 Normals.Append( NewNormals(i) );
648 for (j = 1; j <= Normals.Length(); j++)
650 if (Normals(j).SquareMagnitude() == 0.)
652 Standard_Boolean isFound = Standard_False;
653 for (k = 1; k <= Bset.Length(); k++)
654 if (j == Bset(k).Ind1 || j == Bset(k).Ind2)
656 isFound = Standard_True;
660 Normals(j) = NullVec;
664 return Standard_True;