src/libengrid/geometrytools.h

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