src/SMESHGUI/SMESHGUI_ComputeScalarValue.cxx

   1 using namespace std;
   2 //  File      : SMESHGUI_ComputeScalarValue.cxx
   3 //  Created   : Mon Jun 24 14:06:00 2002
   4 //  Author    : Nicolas REJNERI
   5
   6 //  Project   : SALOME
   7 //  Module    : SMESH
   8 //  Copyright : Open CASCADE 2002
   9 //  $Header$
  10
  11 #include "SMESHGUI_ComputeScalarValue.h"
  12 #include "utilities.h"
  13 #include <math.h>
  14
  15
  16
  17
  18 //=============================================================================
  19 /*!
  20  *
  21  */
  22 //=============================================================================
  23 static double ComputeLength(float* p1, float* p2) {
  24   float a1,a2,a3,b1,b2,b3;
  25   a1 =  p1[0];
  26   a2 =  p1[1];
  27   a3 =  p1[2];
  28   b1 =  p2[0];
  29   b2 =  p2[1];
  30   b3 =  p2[2];
  31   // MESSAGE( a1 << " "<< a2 << " "<< a3 << " " << b1 << " "<< b2 << " "<< b3 );
  32   float X1,Y1,Z1,X2,Z2,Y2;
  33   X1 = b1 - a1;
  34   Y1 = b2 - a2;
  35   Z1 = b3 - a3;
  36   // MESSAGE( X1 << " "<< Y1 << " "<< Z1 );
  37   float e1;
  38   e1 = sqrt( X1*X1 + Y1*Y1 + Z1*Z1 ) ;
  39   // MESSAGE( "Length = " << e1 );
  40   return e1;
  41 }
  42
  43 //=============================================================================
  44 /*!
  45  *
  46  */
  47 //=============================================================================
  48 double SMESHGUI_ComputeScalarValue::LengthEdges(vtkCell* theCell) {
  49   int num_points = theCell->GetNumberOfPoints ();
  50   vtkPoints* points = theCell->GetPoints();
  51   if (num_points != 2 ) return 0;
  52   float* pnt1 = points->GetPoint(0);
  53   float* pnt2 = points->GetPoint(1);
  54   return ComputeLength(pnt1,pnt2);
  55 };
  56
  57 //=============================================================================
  58 /*!
  59  *
  60  */
  61 //=============================================================================
  62 static double ComputeAreaOfTriangle(float* p1, float* p2, float* p3) {
  63   double a1,a2,a3,b1,b2,b3,c1,c2,c3;
  64   a1 =  p1[0];
  65   a2 =  p1[1];
  66   a3 =  p1[2];
  67   b1 =  p2[0];
  68   b2 =  p2[1];
  69   b3 =  p2[2];
  70   c1 =  p3[0];
  71   c2 =  p3[1];
  72   c3 =  p3[2];
  73
  74   float e1, e2, e3;
  75   e1 = sqrt( (a1-b1)*(a1-b1) + (a2-b2)*(a2-b2) + (a3-b3)*(a3-b3) ) ;
  76   e2 = sqrt( (b1-c1)*(b1-c1) + (b2-c2)*(b2-c2) + (b3-c3)*(b3-c3) ) ;
  77   e3 = sqrt( (c1-a1)*(c1-a1) + (c2-a2)*(c2-a2) + (c3-a3)*(c3-a3) ) ;
  78
  79   // MESSAGE( "e = " << e1 << " " << e2 <<"  " << e3 );
  80   float s = (e1+e2+e3)/2;
  81   double area = sqrt(s*(s-e1)*(s-e2)*(s-e3));
  82   // MESSAGE( "area = " << area );
  83   return area;
  84 }
  85
  86 //=============================================================================
  87 /*!
  88  *
  89  */
  90 //=============================================================================
  91 double SMESHGUI_ComputeScalarValue::AreaElements(vtkCell* theCell){
  92   // MESSAGE ( " SMESHGUI_ComputeScalarValue::AreaElements " )
  93   int num_points = theCell->GetNumberOfPoints ();
  94   vtkPoints* points = theCell->GetPoints();
  95   // MESSAGE( "num_points =  "<< num_points );
  96   for (int j = 0; j < theCell->GetNumberOfPoints (); j++) {
  97     float* pnt = points->GetPoint(j);
  98     // MESSAGE( pnt[0] << " " << pnt[1] << " " <<  pnt[2] );
  99   }
 100   if (num_points < 3 ) return 0;
 101   if (num_points == 3) {
 102     float* p1 = points->GetPoint(0);
 103     float* p2 = points->GetPoint(1);
 104     float* p3 = points->GetPoint(2);
 105     double area = ComputeAreaOfTriangle(p1,p2,p3);
 106     return area;
 107   } else if (num_points == 4) {
 108     float* p1 = points->GetPoint(0);
 109     float* p2 = points->GetPoint(1);
 110     float* p3 = points->GetPoint(2);
 111     float* p4 = points->GetPoint(3);
 112     double area1 = ComputeAreaOfTriangle(p1,p2,p3);
 113     double area2 = ComputeAreaOfTriangle(p3,p4,p1);
 114     return area1+area2;
 115   }
 116 };
 117
 118 //=============================================================================
 119 /*!
 120  *
 121  */
 122 //=============================================================================
 123 double  SMESHGUI_ComputeScalarValue::Taper(vtkCell* theCell){
 124   int num_points = theCell->GetNumberOfPoints ();
 125   vtkPoints* points = theCell->GetPoints();
 126   if (num_points != 4 ) return 0;
 127   float* p1 = points->GetPoint(0);
 128   float* p2 = points->GetPoint(1);
 129   float* p3 = points->GetPoint(2);
 130   float* p4 = points->GetPoint(3);
 131   double A1 = ComputeAreaOfTriangle(p4,p1,p2);
 132   double A2 = ComputeAreaOfTriangle(p3,p1,p2);
 133   double A3 = ComputeAreaOfTriangle(p2,p3,p4);
 134   double A4 = ComputeAreaOfTriangle(p3,p4,p1);
 135   double JA = 0.25 * (A1 + A2 + A3 + A4);
 136   double taper = fabs(A1/(JA-1));
 137   if (fabs(A2/(JA-1)) > taper) taper = fabs(A2/(JA-1));
 138   if (fabs(A3/(JA-1)) > taper) taper = fabs(A3/(JA-1));
 139   if (fabs(A4/(JA-1)) > taper) taper = fabs(A4/(JA-1));
 140   // MESSAGE( "Taper = " << taper);
 141   return taper;
 142 };
 143
 144 //=============================================================================
 145 /*!
 146  *
 147  */
 148 //=============================================================================
 149 double  SMESHGUI_ComputeScalarValue::AspectRatio(vtkCell* theCell) {
 150   int num_points = theCell->GetNumberOfPoints ();
 151   vtkPoints* points = theCell->GetPoints();
 152   if (num_points < 3 ) return 0;
 153   if (num_points == 3) {
 154     float a1,a2,a3,b1,b2,b3,c1,c2,c3;
 155     float* pnt = points->GetPoint(0);
 156     a1 =  pnt[0];
 157     a2 =  pnt[1];
 158     a3 =  pnt[2];
 159     pnt = points->GetPoint(1);
 160     b1 =  pnt[0];
 161     b2 =  pnt[1];
 162     b3 =  pnt[2];
 163     pnt = points->GetPoint(2);
 164     c1 =  pnt[0];
 165     c2 =  pnt[1];
 166     c3 =  pnt[2];
 167
 168     float e1, e2, e3;
 169     e1 = sqrt( (a1-b1)*(a1-b1) + (a2-b2)*(a2-b2) + (a3-b3)*(a3-b3) ) ;
 170     e2 = sqrt( (b1-c1)*(b1-c1) + (b2-c2)*(b2-c2) + (b3-c3)*(b3-c3) ) ;
 171     e3 = sqrt( (c1-a1)*(c1-a1) + (c2-a2)*(c2-a2) + (c3-a3)*(c3-a3) ) ;
 172
 173     float amax,p,s;
 174
 175     amax=e1;
 176     if (e2>amax) amax=e2;
 177     if (e3>amax) amax=e3;
 178
 179     p=(e1+e2+e3)/2;
 180     s=AreaElements(theCell);
 181
 182     double aspectRatio=amax*p*sqrt(double(3))/(s*6);
 183     // MESSAGE( "aspectRatio = " << aspectRatio );
 184     return(aspectRatio);
 185   }
 186   else if (num_points == 4) {
 187     float a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3;
 188     float* pnt = points->GetPoint(0);
 189     a1 =  pnt[0];
 190     a2 =  pnt[1];
 191     a3 =  pnt[2];
 192     pnt = points->GetPoint(1);
 193     b1 =  pnt[0];
 194     b2 =  pnt[1];
 195     b3 =  pnt[2];
 196     pnt = points->GetPoint(2);
 197     c1 =  pnt[0];
 198     c2 =  pnt[1];
 199     c3 =  pnt[2];
 200     pnt = points->GetPoint(3);
 201     d1 =  pnt[0];
 202     d2 =  pnt[1];
 203     d3 =  pnt[2];
 204
 205     float e1, e2, e3, e4;
 206     float len_min, len_max;
 207     e1 = sqrt( (a1-b1)*(a1-b1) + (a2-b2)*(a2-b2) + (a3-b3)*(a3-b3) ) ;
 208     e2 = sqrt( (b1-c1)*(b1-c1) + (b2-c2)*(b2-c2) + (b3-c3)*(b3-c3) ) ;
 209     e3 = sqrt( (c1-d1)*(c1-d1) + (c2-d2)*(c2-d2) + (c3-d3)*(c3-d3) ) ;
 210     e4 = sqrt( (d1-a1)*(d1-a1) + (d2-a2)*(d2-a2) + (d3-a3)*(d3-a3) ) ;
 211
 212     len_min = e1; len_max = e1;
 213
 214     if (e2 >len_max ) len_max = e2;
 215     if (e3 >len_max ) len_max = e3;
 216     if (e4 >len_max ) len_max = e4;
 217     if (e2 <len_min ) len_min = e2;
 218     if (e3 <len_min ) len_min = e3;
 219     if (e4 <len_min ) len_min = e4;
 220
 221     return (len_max/len_min);
 222   }
 223 };
 224
 225 //=============================================================================
 226 /*!
 227  *
 228  */
 229 //=============================================================================
 230 static double ComputeAngle(float* p1, float* p2, float* p3) {
 231   const double pi=4*atan(double(1));
 232   float a1,a2,a3,b1,b2,b3,c1,c2,c3;
 233   a1 =  p1[0];
 234   a2 =  p1[1];
 235   a3 =  p1[2];
 236   b1 =  p2[0];
 237   b2 =  p2[1];
 238   b3 =  p2[2];
 239   c1 =  p3[0];
 240   c2 =  p3[1];
 241   c3 =  p3[2];
 242   float X1,Y1,Z1,X2,Z2,Y2;
 243   X1 = b1 - a1;
 244   X2 = c1 - b1;
 245   Y1 = b2 - a2;
 246   Y2 = c2 - b2;
 247   Z1 = b3 - a3;
 248   Z2 = c3 - b3;
 249
 250
 251   float e1, e2, e3;
 252   e1 = sqrt( X1*X1 + Y1*Y1 + Z1*Z1 ) ;
 253   e2 = sqrt( X2*X2 + Y2*Y2 + Z2*Z2 ) ;
 254   double dot=(X1*(X2)+Y1*(Y2)+Z1*(Z2));
 255   //  MESSAGE( pi );
 256   //  MESSAGE( dot/(e1*e2) );
 257   double cosinus = dot/(e1*e2);
 258   cosinus = fabs(cosinus);
 259   return 180*acos (cosinus)/pi;
 260 }
 261
 262 //=============================================================================
 263 /*!
 264  *
 265  */
 266 //=============================================================================
 267 double SMESHGUI_ComputeScalarValue::MinimumAngle(vtkCell* theCell) {
 268   int num_points = theCell->GetNumberOfPoints ();
 269   vtkPoints* points = theCell->GetPoints();
 270   if (num_points < 3 ) return 0;
 271   float* pnt1 = points->GetPoint(0);
 272   float* pnt2 = points->GetPoint(1);
 273   float* pnt3 = points->GetPoint(2);
 274   if (num_points == 3) {
 275     double a1,a2,a3,amin;
 276     a1=fabs(ComputeAngle(pnt1,pnt2,pnt3));
 277     amin=a1;
 278     a2=fabs(ComputeAngle(pnt2,pnt3,pnt1));
 279     if (a2<amin) amin=a2;
 280     a3=fabs(ComputeAngle(pnt3,pnt1,pnt2));
 281     if (a3<amin) amin=a3;
 282     // MESSAGE( "Minimal angle " << amin );
 283     return amin;
 284   }
 285   else if (num_points == 4) {
 286     float* pnt4 = points->GetPoint(3);
 287     double a1,a2,a3,a4,amin;
 288     a1=fabs(ComputeAngle(pnt1,pnt2,pnt3));
 289     amin=a1;
 290     a2=fabs(ComputeAngle(pnt2,pnt3,pnt4));
 291     if (a2<amin) amin=a2;
 292     a3=fabs(ComputeAngle(pnt3,pnt4,pnt1));
 293     if (a3<amin) amin=a3;
 294     a4=fabs(ComputeAngle(pnt4,pnt1,pnt2));
 295     if (a4<amin) amin=a4;
 296
 297     // MESSAGE( "Minimal angle " << amin );
 298     return amin;
 299   }
 300 };
 301
 302 //=============================================================================
 303 /*!
 304  *
 305  */
 306 //=============================================================================
 307 double  SMESHGUI_ComputeScalarValue::Skew(vtkCell* theCell) {
 308   int num_points = theCell->GetNumberOfPoints ();
 309   vtkPoints* points = theCell->GetPoints();
 310   if (num_points < 3 ) return 0;
 311   //triangle case
 312   if (num_points == 3) {
 313     float* pnt1 = points->GetPoint(0);
 314     float* pnt2 = points->GetPoint(1);
 315     float* pnt3 = points->GetPoint(2);
 316     double a1,a2,a3,amax;
 317     a1=fabs(60 - fabs(ComputeAngle(pnt1,pnt2,pnt3)));
 318     amax=a1;
 319     a2=fabs(60 - fabs(ComputeAngle(pnt2,pnt3,pnt1)));
 320     if (a2>amax) amax=a2;
 321     a3=fabs(60 - fabs(ComputeAngle(pnt3,pnt1,pnt2)));
 322     if (a3>amax) amax=a3;
 323     // MESSAGE( "Skew = " << amax );
 324     return amax;
 325   }
 326   //quadrangle case
 327   else if (num_points == 4) {
 328     float* pnt1 = points->GetPoint(0);
 329     float* pnt2 = points->GetPoint(1);
 330     float* pnt3 = points->GetPoint(2);
 331     float* pnt4 = points->GetPoint(3);
 332
 333     double a1,a2,a3,a4,amax;
 334     a1=fabs(90 - fabs(ComputeAngle(pnt1,pnt2,pnt3)));
 335     amax=a1;
 336     a2=fabs(90 - fabs(ComputeAngle(pnt2,pnt3,pnt4)));
 337     if (a2>amax) amax=a2;
 338     a3=fabs(90 - fabs(ComputeAngle(pnt3,pnt4,pnt1)));
 339     if (a3>amax) amax=a3;
 340     a4=fabs(90 - fabs(ComputeAngle(pnt4,pnt1,pnt2)));
 341     if (a4>amax) amax=a4;
 342     // MESSAGE( "Skew = " << amax );
 343     return amax;
 344   }
 345 };
 346
 347 //=============================================================================
 348 /*!
 349  *
 350  */
 351 //=============================================================================
 352 static double ComputeA(float* p1, float* p2, float* p3, float* G) {
 353   double e1 = sqrt(pow(p2[0]-p1[0], 2)+pow(p2[1]-p1[1], 2)+pow(p2[2]-p1[2], 2));
 354   double e2 = sqrt(pow(p3[0]-p2[0], 2)+pow(p3[1]-p2[1], 2)+pow(p3[2]-p2[2], 2));
 355   double l;
 356   if (e1 < e2) l = 0.5*e1;
 357   else l = 0.5*e2;
 358   float GI[3], GJ[3], N[3];;
 359   GI[0] = (p2[0]-p1[0])/2-G[0];
 360   GI[1] = (p2[1]-p1[1])/2-G[1];
 361   GI[2] = (p2[2]-p1[2])/2-G[2];
 362
 363   GJ[0] = (p3[0]-p2[0])/2-G[0];
 364   GJ[1] = (p3[1]-p2[1])/2-G[1];
 365   GJ[2] = (p3[2]-p2[2])/2-G[2];
 366
 367   N[0] = GI[1]*GJ[2] - GI[2]*GJ[1];
 368   N[1] = GI[2]*GJ[0] - GI[0]*GJ[2];
 369   N[2] = GI[0]*GJ[1] - GI[1]*GJ[0];
 370
 371   double H;
 372   float T[3];
 373   T[0] = (p1[0]-G[0])*N[0];
 374   T[1] = (p1[1]-G[1])*N[1];
 375   T[2] = (p1[2]-G[2])*N[2];
 376
 377   H = sqrt(pow(T[0],2)+pow(T[1],2)+pow(T[2],2))/sqrt(pow(N[0],2)+pow(N[1],2)+pow(N[2],2));
 378   double A = H/l;
 379   return A;
 380 }
 381
 382 //=============================================================================
 383 /*!
 384  *
 385  */
 386 //=============================================================================
 387 double  SMESHGUI_ComputeScalarValue::Warp(vtkCell* theCell) {
 388   int num_points = theCell->GetNumberOfPoints ();
 389   vtkPoints* points = theCell->GetPoints();
 390   if (num_points != 4 ) return 0;
 391   float* p1 = points->GetPoint(0);
 392   float* p2 = points->GetPoint(1);
 393   float* p3 = points->GetPoint(2);
 394   float* p4 = points->GetPoint(3);
 395   double G1, G2, G3;
 396   float G[3];
 397   G[0] = (p1[0]+p2[0]+p3[0]+p4[0])/4;
 398   G[1] = (p1[1]+p2[1]+p3[1]+p4[1])/4;
 399   G[2] = (p1[2]+p2[2]+p3[2]+p4[2])/4;
 400   double amax = ComputeA(p1, p2, p3, G);
 401   double nextA = ComputeA(p2, p3, p4, G);
 402   if (nextA > amax) amax = nextA;
 403   nextA = ComputeA(p3, p4, p1, G);
 404   if (nextA > amax) amax = nextA;
 405   nextA = ComputeA(p4, p1, p2, G);
 406   if (nextA > amax) amax = nextA;
 407   // MESSAGE( "Warp = " << amax );
 408   return amax;
 409 }