src/IntAna2d/IntAna2d_Outils.cxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 //============================================================================
  16 //======================================================= IntAna2d_Outils.cxx
  17 //============================================================================
  18 #include <IntAna2d_Outils.hxx>
  19 #include <math_DirectPolynomialRoots.hxx>
  20
  21 MyDirectPolynomialRoots::MyDirectPolynomialRoots(const Standard_Real A4,
  22                                                  const Standard_Real A3,
  23                                                  const Standard_Real A2,
  24                                                  const Standard_Real A1,
  25                                                  const Standard_Real A0) {
  26   //-- std::cout<<" IntAna2d : A4..A0 "<<A4<<" "<<A3<<" "<<A2<<" "<<A1<<" "<<A0<<" "<<std::endl;
  27   nbsol = 0;
  28   same = Standard_False;
  29 //  Modified by Sergey KHROMOV - Thu Oct 24 13:10:14 2002 Begin
  30   Standard_Real anAA[5];
  31
  32   anAA[0] = Abs(A0);
  33   anAA[1] = Abs(A1);
  34   anAA[2] = Abs(A2);
  35   anAA[3] = Abs(A3);
  36   anAA[4] = Abs(A4);
  37
  38 //   if((Abs(A4)+Abs(A3)+Abs(A2)+Abs(A1)+Abs(A0))<Epsilon(10000.0))  {
  39   if((anAA[0]+anAA[1]+anAA[2]+anAA[3]+anAA[4])<Epsilon(10000.0))  {
  40 //  Modified by Sergey KHROMOV - Thu Oct 24 13:10:15 2002 End
  41     same = Standard_True;
  42     return;
  43   }
  44   Standard_Integer i,j,nbp;
  45   for(i=0;i<16;i++) val[i]=RealLast();
  46
  47   Standard_Real tol = Epsilon(100.0);
  48   math_DirectPolynomialRoots MATH_A43210(A4,A3,A2,A1,A0);
  49   Standard_Boolean PbPossible = Standard_False;
  50   Standard_Integer NbsolPolyComplet = 0;
  51   if(MATH_A43210.IsDone()) {
  52     nbp = MATH_A43210.NbSolutions();
  53     NbsolPolyComplet = nbp;
  54     for(i=1;i<=nbp;i++) {
  55       Standard_Real x = MATH_A43210.Value(i);
  56       //-- std::cout<<" IntAna2d : x Pol Complet :"<<x<<std::endl;
  57       val[nbsol] = A0 + x*(A1+x*(A2+x*(A3+x*A4)));
  58       sol[nbsol] = x;
  59       if(val[nbsol]> tol  ||  val[nbsol]<-tol) {
  60         PbPossible = Standard_True;
  61       }
  62       nbsol++;
  63     }
  64     if(nbp & 1)
  65       PbPossible = Standard_True;
  66   }
  67   else {
  68     PbPossible = Standard_True;
  69   }
  70   //-- On recherche le plus petit coeff entre A4 et A0
  71   if(PbPossible) {
  72 //  Modified by Sergey KHROMOV - Thu Oct 24 12:45:35 2002 Begin
  73     Standard_Real anAMin = RealLast();
  74     Standard_Real anAMax = -1;
  75     Standard_Real anEps  = RealEpsilon();
  76
  77     for (i = 0; i < 5; i++) {
  78       anAMin = Min(anAMin, Max(anAA[i], anEps));
  79       anAMax = Max(anAMax, Max(anAA[i], anEps));
  80     }
  81
  82     anEps = Min(1.e-4, Epsilon(1000.*anAMax/anAMin));
  83 //  Modified by Sergey KHROMOV - Thu Oct 24 15:46:24 2002 End
  84     math_DirectPolynomialRoots MATH_A4321(A4,A3,A2,A1);
  85     if(MATH_A4321.IsDone()) {
  86       nbp = MATH_A4321.NbSolutions();
  87       //-- On Ajoute les valeurs au tableau
  88       for(i=1;i<=nbp;i++) {
  89         Standard_Real x = MATH_A4321.Value(i);
  90         Standard_Boolean Add = Standard_True;
  91         for(j=0;j<nbsol;j++) {
  92           Standard_Real t = sol[j]-x;
  93 //  Modified by Sergey KHROMOV - Thu Oct 24 12:04:26 2002 Begin
  94 //        if(Abs(t)<tol) {
  95           if(Abs(t) < anEps) {
  96 //  Modified by Sergey KHROMOV - Thu Oct 24 12:04:47 2002 End
  97             Add = Standard_False;
  98           }
  99         }
 100         if(Add) {
 101           val[nbsol] =  A0 + x*(A1+x*(A2+x*(A3+x*A4)));
 102           sol[nbsol] = x;
 103           nbsol++;
 104         }
 105       }
 106     }
 107     math_DirectPolynomialRoots MATH_A3210(A3,A2,A1,A0);
 108     if(MATH_A3210.IsDone()) {
 109       nbp = MATH_A3210.NbSolutions();
 110       //-- On Ajoute les valeurs au tableau
 111       for(i=1;i<=nbp;i++) {
 112         Standard_Real x = MATH_A3210.Value(i);
 113         Standard_Boolean Add = Standard_True;
 114         for(j=0;j<nbsol;j++) {
 115           Standard_Real t = sol[j]-x;
 116 //  Modified by Sergey KHROMOV - Thu Oct 24 12:06:01 2002 Begin
 117 //        if(Abs(t)<tol) {
 118           if(Abs(t) < anEps) {
 119 //  Modified by Sergey KHROMOV - Thu Oct 24 12:06:04 2002 End
 120             Add = Standard_False;
 121           }
 122         }
 123         if(Add) {
 124           val[nbsol] =  A0 + x*(A1+x*(A2+x*(A3+x*A4)));
 125           sol[nbsol] = x;
 126           nbsol++;
 127         }
 128       }
 129     }
 130     math_DirectPolynomialRoots MATH_A210(A3,A2,A1);
 131     if(MATH_A210.IsDone()) {
 132       nbp = MATH_A210.NbSolutions();
 133       //-- On Ajoute les valeurs au tableau
 134       for(i=1;i<=nbp;i++) {
 135         Standard_Real x = MATH_A210.Value(i);
 136         Standard_Boolean Add = Standard_True;
 137         for(j=0;j<nbsol;j++) {
 138           Standard_Real t = sol[j]-x;
 139 //  Modified by Sergey KHROMOV - Thu Oct 24 12:07:04 2002 Begin
 140 //        if(Abs(t)<tol) {
 141           if(Abs(t) < anEps) {
 142 //  Modified by Sergey KHROMOV - Thu Oct 24 12:07:06 2002 End
 143             Add = Standard_False;
 144           }
 145         }
 146         if(Add) {
 147           val[nbsol] =  A0 + x*(A1+x*(A2+x*(A3+x*A4)));
 148           sol[nbsol] = x;
 149           nbsol++;
 150         }
 151       }
 152     }
 153     //------------------------------------------------------------
 154     //-- On trie les valeurs par ordre decroissant de val
 155     //-- for(i=0;i<nbsol;i++) {
 156     //--  std::cout<<" IntAna2d Sol,Val"<<sol[i]<<"  "<<val[i]<<std::endl;
 157     //-- }
 158     Standard_Boolean TriOK = Standard_False;
 159     do {
 160       TriOK = Standard_True;
 161       for(i=1; i<nbsol;i++) {
 162         if(Abs(val[i])<Abs(val[i-1])) {
 163           Standard_Real t;
 164           t        = val[i];
 165           val[i]   = val[i-1];
 166           val[i-1] = t;
 167           t        = sol[i];
 168           sol[i]   = sol[i-1];
 169           sol[i-1] = t;
 170           TriOK = Standard_False;
 171         }
 172       }
 173     }
 174     while(!TriOK);
 175     //-----------------------------------------------------------
 176     //-- On garde les premieres valeurs
 177     //-- Au moins autant que le polynome Complet
 178     //--
 179     for(nbsol=0; nbsol<NbsolPolyComplet || Abs(val[nbsol])<Epsilon(10000.0); nbsol++) ;
 180     //-- std::cout<<" IntAna2d : nbsol:"<<nbsol<<std::endl;
 181   }
 182   if(nbsol==0) {
 183     nbsol=-1;
 184   }
 185   if(nbsol>4) {
 186     same=1;
 187     nbsol=0;
 188   }
 189 }
 190
 191
 192 MyDirectPolynomialRoots::MyDirectPolynomialRoots(const Standard_Real A2,
 193                                                  const Standard_Real A1,
 194                                                  const Standard_Real A0) {
 195   //-- std::cout<<" IntAna2d : A2..A0 "<<A2<<" "<<A1<<" "<<A0<<" "<<std::endl;
 196   nbsol=0;
 197   if((Abs(A2)+Abs(A1)+Abs(A0))<Epsilon(10000.0))  {
 198     same = Standard_True;
 199     return;
 200   }
 201   math_DirectPolynomialRoots MATH_A210(A2,A1,A0);
 202   same = Standard_False;
 203   if(MATH_A210.IsDone()) {
 204     for(Standard_Integer i=1;i<=MATH_A210.NbSolutions(); i++) {
 205       Standard_Real x = MATH_A210.Value(i);
 206       val[nbsol] = A0 + x*(A1+x*A2);
 207       sol[nbsol] = x;
 208       //-- std::cout<<" IntAna2d : x Pol Complet :"<<x<<"  Val:"<<val[nbsol]<<std::endl;
 209       nbsol++;
 210     }
 211   }
 212   else {
 213     nbsol = -1;
 214   }
 215 }
 216
 217 Standard_Boolean Points_Confondus(const Standard_Real x1,const Standard_Real y1
 218                                   ,const Standard_Real x2,const Standard_Real y2) {
 219   if(Abs(x1-x2)<Epsilon(x1)) {
 220     if(Abs(y1-y2)<Epsilon(y1)) {
 221       return(Standard_True);
 222     }
 223   }
 224   return(Standard_False);
 225 }
 226
 227 //-----------------------------------------------------------------------------
 228 //--- Les points confondus sont supprimes
 229 //--- Le nombre de points est mis a jour
 230
 231 void Traitement_Points_Confondus(Standard_Integer& nb_pts,
 232                                  IntAna2d_IntPoint* pts) {
 233   Standard_Integer i,j;
 234   for(i=nb_pts;i>1;i--) {
 235     Standard_Boolean Non_Egalite=Standard_True;
 236     for(j=i-1;(j>0) && Non_Egalite;j--) {
 237       //                                        <--- Deja Teste --->
 238       //             | 1  |2  |  | J |  |I-1| I |I+1|          |NPTS|
 239       //             | 1  |2  |  | J |  |I-1|XXX|I+1|          |NPTS|
 240       //             | 1  |2  |  | J |  |I-1|I+1|I+2|     |NPTS|
 241       if(Points_Confondus(pts[i-1].Value().X(),
 242                           pts[i-1].Value().Y(),
 243                           pts[j-1].Value().X(),
 244                           pts[j-1].Value().Y())) {
 245         Standard_Integer k;
 246         Non_Egalite=Standard_False;
 247         for(k=i;k<nb_pts;k++) {
 248           Standard_Real Xk,Yk,Uk;
 249           Xk=pts[k].Value().X();
 250           Yk=pts[k].Value().Y();
 251           Uk=pts[k].ParamOnFirst();
 252           pts[k-1].SetValue(Xk,Yk,Uk);
 253         }
 254         nb_pts--;
 255       }
 256     }
 257   }
 258 }
 259
 260 //-----------------------------------------------------------------------------
 261 void Coord_Ancien_Repere(Standard_Real& x1,Standard_Real& y1,const gp_Ax2d Dir1) {
 262   Standard_Real t11,t12,t21,t22,t13,t23;
 263   Standard_Real x0,y0;
 264
 265   // x1 et y1 Sont les Coordonnees dans le repere lie a Dir1
 266   // On Renvoie ces Coordonnees dans le repere "absolu"
 267
 268   Dir1.Direction().Coord(t11,t21);
 269   Dir1.Location().Coord(t13,t23);
 270
 271   t22=t11;
 272   t12=-t21;
 273
 274   x0= t11*x1 + t12*y1 + t13;
 275   y0= t21*x1 + t22*y1 + t23;
 276
 277   x1=x0;
 278   y1=y0;
 279 }
 280
 281
 282
 283 #if 0
 284
 285 //-- A Placer dans les ressources de la classe Conic   ??
 286 //-----------------------------------------------------------------------------
 287 //--- Calcul des Coefficients A,..F dans le repere lie a  Dir1
 288 //--- A Partir des Coefficients dans le repere "Absolu"
 289
 290 void Coeff_Nouveau_Repere(Standard_Real& A,Standard_Real& B,Standard_Real& C
 291                           ,Standard_Real& D,Standard_Real& E,Standard_Real& F
 292                           ,const gp_Ax2d Dir1)  {
 293   Standard_Real t11,t12,t13;                  // x = t11 X + t12 Y + t13
 294   Standard_Real t21,t22,t23;                  // y = t21 X + t22 Y + t23
 295   Standard_Real A1,B1,C1,D1,E1,F1;
 296
 297   // On a P0(x,y)=A x x + B y y + ... + F =0    (x et y ds le repere "Absolu")
 298   // et on cherche P1(X(x,y),Y(x,y))=P0(x,y)
 299   // Avec P1(X,Y)= A1 X X + B1 Y Y + 2 C1 X Y + 2 D1 X + 2 E1 Y + F1
 300   //             = A  x x + B  y y + 2 C  x y + 2 D  x + 2 E  y + f
 301
 302   Dir1.Direction().Coord(t11,t21);
 303   Dir1.Location().Coord(t13,t23);
 304
 305   t22=t11;
 306   t12=-t21;
 307
 308   A1=(t11*(A*t11 + 2*C*t21) + B*t21*t21);
 309   B1=(t12*(A*t12 + 2*C*t22) + B*t22*t22);
 310   C1=(t12*(A*t11 + C*t21) + t22*(C*t11 + B*t21));
 311   D1=(t11*(D + A*t13) + t21*(E + C*t13) + t23*(C*t11 + B*t21));
 312   E1=(t12*(D + A*t13) + t22*(E + C*t13) + t23*(C*t12 + B*t22));
 313   F1=F + t13*(2.0*D + A*t13) + t23*(2.0*E + 2.0*C*t13 + B*t23);
 314
 315   A=A1; B=B1; C=C1; D=D1; E=E1; F=F1;
 316 }
 317 #endif
 318
 319
 320
 321
 322
 323
 324
 325
 326
 327