src/libengrid/geometrytools.h

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2012 enGits GmbH                                     +
   7 // +                                                                      +
   8 // + enGrid is free software: you can redistribute it and/or modify       +
   9 // + it under the terms of the GNU General Public License as published by +
  10 // + the Free Software Foundation, either version 3 of the License, or    +
  11 // + (at your option) any later version.                                  +
  12 // +                                                                      +
  13 // + enGrid is distributed in the hope that it will be useful,            +
  14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  16 // + GNU General Public License for more details.                         +
  17 // +                                                                      +
  18 // + You should have received a copy of the GNU General Public License    +
  19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  20 // +                                                                      +
  21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  22 //
  23 #ifndef geometrytools_H
  24 #define geometrytools_H
  25
  26 #include "math/mathvector.h"
  27 #include "math/smallsquarematrix.h"
  28
  29 #include <QVector>
  30
  31 #include <vtkUnstructuredGrid.h>
  32
  33 namespace GeometryTools
  34 {
  35
  36
  37 double rad2deg( double rad ); ///< Converts radians to degrees
  38 double deg2rad( double deg ); ///< Converts degrees to radians
  39
  40 void rotate(vec3_t g1, vec3_t g2, vec3_t g3, vec3_t &b, double theta);
  41
  42 /** Rotates vector v around vector axis by an angle theta */
  43 vec3_t rotate(vec3_t v, vec3_t axis, double theta);
  44
  45 /** Returns a normalized vector orthogonal to v */
  46 vec3_t orthogonalVector(vec3_t v);
  47
  48 /** Calculates the intersection between a line and a plane.
  49  * @param x_straight A point of the line.
  50  * @param v_straight Direction of the line.
  51  * @param x_plane A point of the plane
  52  * @param n_plane Normal of the plane
  53  */
  54 double intersection(vec3_t x_straight, vec3_t v_straight,
  55                     vec3_t x_plane, vec3_t n_plane);
  56
  57 /** Calculates the intersection between a line and a plane.
  58  * @param x_straight A point of the line.
  59  * @param v_straight Direction of the line.
  60  * @param x_plane A point of the plane
  61  * @param u_plane A vector of the plane
  62  * @param v_plane Another vector of the plane not colinear with u_plane
  63  */
  64 double intersection(vec3_t x_straight, vec3_t v_straight,
  65                     vec3_t x_plane, vec3_t u_plane, vec3_t v_plane);
  66
  67 /**
  68  * Calculates the intersection of a segment [x1,x2] with a triangle (a,b,c)
  69  * Note: (xi,ri) will always be set to the intersection of the line (x1,x2) with the plane (a,b,c) even if the segment does not intersect the triangle!
  70  * @param a Input: Triangle point 1
  71  * @param b Input: Triangle point 2
  72  * @param c Input: Triangle point 3
  73  * @param x1 Input: Segment point 1
  74  * @param x2 Input: Segment point 2
  75  * @param xi Output: 3D Global coordinates of the intersection point
  76  * @param ri Output: 3D local triangle coordinates of the intersection point
  77  * @param tol Input: Relative tolerance: There can only be an intersection if:
  78  * 0-tol*(x1-x2).abs()<=k<=1+tol*(x1-x2).abs()
  79  * 0-tol<=ri[0]<=1+tol
  80  * 0-tol<=ri[1]<=1+tol
  81  * @return true if an intersection point was found, else false.
  82 */
  83 bool intersectEdgeAndTriangle(const vec3_t& a, const vec3_t& b, const vec3_t& c,
  84                               const vec3_t& x1, const vec3_t& x2, vec3_t& xi, vec3_t& ri, double tol = 1e-4);
  85
  86 bool isInsideTriangle(vec2_t t_M, double tol = 1e-4);
  87
  88 /** Calculates the intersection point M between the lines (r1,u1) and (r2,u2).
  89  * @param k1 Returned by reference. Verifies M = r1+k1*u1
  90  * @param k2 Returned by reference. Verifies M = r2+k2*u2
  91  * @param r1 point of line 1
  92  * @param r2 point of line 2
  93  * @param u1 direction of line 1
  94  * @param u2 direction of line 2
  95  * @return true if an intersection has been found, else false.
  96  */
  97 bool intersection (double &k1, double &k2, vec2_t r1, vec2_t u1, vec2_t r2, vec2_t u2);
  98
  99 void sliceTriangle(const vector<vec3_t> &Tin, vec3_t x, vec3_t n, vector<vector<vec3_t> > &Tout);
 100
 101 /** Returns the volume of a tetrahedron.
 102  * V= v1*(v2^v3) with vi=xi-x0
 103  * If neg = false and V<0, it will return V=-1e99, else it returns V.
 104  * @param x0 point 0 of the tetrahedron
 105  * @param x1 point 1 of the tetrahedron
 106  * @param x2 point 2 of the tetrahedron
 107  * @param x3 point 3 of the tetrahedron
 108  * @param neg If neg = false and V<0, it will return V=-1e99, else it returns V.
 109  * @return volume of the tetrahedron
 110  */
 111 double tetraVol(const vec3_t& x0, const vec3_t& x1, const vec3_t& x2, const vec3_t& x3, bool neg = false);
 112
 113 double pyraVol(vec3_t x1, vec3_t x2, vec3_t x3, vec3_t x4, vec3_t x5, bool neg = false);
 114
 115 double prismVol(vec3_t x1, vec3_t x2, vec3_t x3, vec3_t x4, vec3_t x5, vec3_t x6, bool neg = false);
 116
 117 double hexaVol(vec3_t x1, vec3_t x2, vec3_t x3, vec3_t x4, vec3_t x5, vec3_t x6, vec3_t x7, vec3_t x8,  bool neg = false);
 118
 119 double triArea(vec3_t x1, vec3_t x2, vec3_t x3);
 120
 121 double quadArea(vec3_t x1, vec3_t x2, vec3_t x3, vec3_t x4);
 122
 123 vec3_t triNormal(vec3_t x0, vec3_t x1, vec3_t x2);///< Returns the normal of the surface defined by x0,x1,x2
 124
 125 vec3_t quadNormal(vec3_t x0, vec3_t x1, vec3_t x2, vec3_t x3);///< Returns the normal of the surface defined by x0,x1,x2,x3
 126
 127 vec3_t triNormal(vtkUnstructuredGrid *grid, vtkIdType p1, vtkIdType p2, vtkIdType p3);
 128
 129 vec3_t quadNormal(vtkUnstructuredGrid *grid, vtkIdType p1, vtkIdType p2, vtkIdType p3, vtkIdType p4);
 130
 131 vec3_t cellNormal(vtkUnstructuredGrid *grid, vtkIdType i);
 132
 133 /// Returns the area or volume of a cell.
 134 double cellVA(vtkUnstructuredGrid *grid, vtkIdType cellId, bool neg = false);
 135
 136 inline vec2_t turnRight(const vec2_t &v)
 137 {
 138   vec2_t u;
 139   u[0] =  v[1];
 140   u[1] = -v[0];
 141   return u;
 142 }
 143
 144 inline vec2_t turnLeft(const vec2_t &v)
 145 {
 146   vec2_t u;
 147   u[0] = -v[1];
 148   u[1] =  v[0];
 149   return u;
 150 }
 151
 152 //polygon must be numbered clockwise
 153 inline bool IsConvex(vec3_t a,vec3_t b,vec3_t c,vec3_t d)
 154 {
 155   vec3_t u[4];
 156   u[0]=b-a;
 157   u[1]=c-b;
 158   u[2]=d-c;
 159   u[3]=a-d;
 160
 161   for(int i=0;i<4;i++) {
 162     vec3_t n=u[i].cross(u[(i+1)%4]);
 163     if(n[2]>0) return(false);
 164   }
 165   return(true);
 166 }
 167
 168 inline bool IsConvex(vec2_t a_2D,vec2_t b_2D,vec2_t c_2D,vec2_t d_2D)
 169 {
 170   vec3_t a_3D(a_2D[0],a_2D[1]);
 171   vec3_t b_3D(b_2D[0],b_2D[1]);
 172   vec3_t c_3D(c_2D[0],c_2D[1]);
 173   vec3_t d_3D(d_2D[0],d_2D[1]);
 174   return(IsConvex(a_3D,b_3D,c_3D,d_3D));
 175 }
 176
 177 /// return the angle with relation to another 3-vector
 178 double angle(const vec3_t & u, const vec3_t & v);
 179
 180 /** return the deviation p1->p2->p3 (angle(p2-p1,p3-p2)) */
 181 double deviation(vtkUnstructuredGrid *grid, vtkIdType p1, vtkIdType p2, vtkIdType p3);
 182
 183 /** return the angle p1,p2,p3 (angle(p1-p2,p3-p2)) */
 184 double angle(vtkUnstructuredGrid *grid, vtkIdType p1, vtkIdType p2, vtkIdType p3);
 185
 186 /** return the cosine of the angle between the normals of cell1 and cell2 */
 187 double cosAngle(vtkUnstructuredGrid *grid, vtkIdType cell1, vtkIdType cell2);
 188
 189 /** Returns the center of mass of cellId + passes the minimal and maximal center to corner distances by reference */
 190 vec3_t getCenter(vtkUnstructuredGrid *grid, vtkIdType cellId, double& Rmin, double& Rmax);
 191
 192 /** Returns the distance between id_node1 and id_node2 */
 193 double distance(vtkUnstructuredGrid *grid, vtkIdType id_node1, vtkIdType id_node2);
 194
 195 /** Returns the distance squared between id_node1 and id_node2 */
 196 double distance2(vtkUnstructuredGrid *grid, vtkIdType id_node1, vtkIdType id_node2);
 197
 198 /** area of the circumscribed circle of the triangle */
 199 double areaOfCircumscribedCircle(vtkUnstructuredGrid *grid, vtkIdType id_cell);
 200
 201 /** Compute the circumscribed circle of a triangle in 3D coordinates.
 202   * @param a first node of the triangle
 203   * @param b second node of the triangle
 204   * @param c third node of the triangle
 205   * @param x on return this will be the centre of the circumscribed circle
 206   * @param radius on return this will be the radius of the circumscribed circle
 207   */
 208 void computeCircumscribedCircle(vec3_t a, vec3_t b, vec3_t c, vec3_t &x, double &radius);
 209
 210 vec3_t getBarycentricCoordinates(double x, double y);
 211
 212 vec3_t intersectionOnPlane(vec3_t v, vec3_t A, vec3_t nA, vec3_t B, vec3_t nB);
 213
 214 /** Projects vector V onto plane (O,i,j)
 215  * @param V The vector to project
 216  * @param i A vector of the plane
 217  * @param j A vector of the plane
 218  * @return Returns a 2D vector (x,y) so that V = x*i +y*j
 219  */
 220 vec2_t projectVectorOnPlane(vec3_t V,vec3_t i,vec3_t j);
 221
 222 vec3_t projectPointOnEdge(const vec3_t& M,const vec3_t& A, const vec3_t& u);
 223
 224 vec3_t projectPointOnPlane(const vec3_t& M, const vec3_t& A, const vec3_t& N);
 225
 226 void cart2spherical(vec3_t x, double &alpha, double& beta, double& r);
 227
 228 inline void cart2spherical(vec3_t x, const vec3_t& x0, double &alpha, double& beta, double& r)
 229 {
 230   return cart2spherical(x - x0, alpha, beta, r);
 231 }
 232
 233 inline vec3_t spherical2cart(double alpha, double beta, double r)
 234 {
 235   return r*vec3_t(cos(alpha)*cos(beta), sin(alpha)*cos(beta), sin(beta));
 236 }
 237
 238 inline vec3_t spherical2cart(vec3_t x0, double alpha, double beta, double r)
 239 {
 240   return x0 + spherical2cart(alpha, beta, r);
 241 }
 242
 243 template <class C>
 244 void planeFit(const C &pts, vec3_t &x0, vec3_t &n, bool closed_loop = false)
 245 {
 246   x0 = vec3_t(0,0,0);
 247   int N = pts.size();
 248   if (!closed_loop) {
 249     ++N;
 250   }
 251   QVector<vec3_t> x(N);
 252   for (int i = 0; i < pts.size(); ++i) {
 253     x[i] = pts[i];
 254     x0 += x[i];
 255   }
 256   if (!closed_loop) {
 257     x.last() = x.first();
 258   }
 259   x0 *= 1.0/pts.size();
 260   for (int i = 0; i < x.size(); ++i) {
 261     x[i] -= x0;
 262   }
 263   n = vec3_t(0,0,0);
 264   for (int i = 0; i < x.size() - 1; ++i) {
 265     n += x[i].cross(x[i+1]);
 266   }
 267   double a11 = 0, a12 = 0, a13 = 0;
 268   double a21 = 0, a22 = 0, a23 = 0;
 269   double a31 = 0, a32 = 0, a33 = 0;
 270   for (int i = 0; i < x.size() - 1; ++i) {
 271     a11 += x[i][0]*x[i][0];
 272     a12 += x[i][0]*x[i][1];
 273     a13 += x[i][0]*x[i][2];
 274
 275     a21 += x[i][0]*x[i][1];
 276     a22 += x[i][1]*x[i][1];
 277     a23 += x[i][1]*x[i][2];
 278
 279     a31 += x[i][0]*x[i][2];
 280     a32 += x[i][1]*x[i][2];
 281     a33 += x[i][2]*x[i][2];
 282   }
 283   //n[2] = n[2];
 284   //n[1] = -a33*n[2]/(a22*a31 - a21*a32);
 285   //n[0] = -(a12*n[1] + a13*n[2])/a11;
 286   n.normalise();
 287 }
 288
 289 template <class C>
 290 vec3_t polyNormal(const C &pts, bool closed_loop = false)
 291 {
 292   vec3_t x0(0,0,0);
 293   int N = pts.size();
 294   if (!closed_loop) {
 295     ++N;
 296   }
 297   QVector<vec3_t> x(N);
 298   for (int i = 0; i < pts.size(); ++i) {
 299     x[i] = pts[i];
 300     x0 += x[i];
 301   }
 302   if (!closed_loop) {
 303     x.last() = x.first();
 304   }
 305   x0 *= 1.0/pts.size();
 306   for (int i = 0; i < x.size(); ++i) {
 307     x[i] -= x0;
 308   }
 309   vec3_t n(0,0,0);
 310   for (int i = 0; i < x.size() - 1; ++i) {
 311     n += x[i].cross(x[i+1]);
 312   }
 313   return n;
 314 }
 315
 316 }
 317
 318 #endif