simgear/scene/sky/clouds3d/tri.cpp

   1 //------------------------------------------------------------------------------
   2 // File : tri.cpp
   3 //------------------------------------------------------------------------------
   4 // GLVU : Copyright 1997 - 2002
   5 //        The University of North Carolina at Chapel Hill
   6 //------------------------------------------------------------------------------
   7 // Permission to use, copy, modify, distribute and sell this software and its
   8 // documentation for any purpose is hereby granted without fee, provided that
   9 // the above copyright notice appear in all copies and that both that copyright
  10 // notice and this permission notice appear in supporting documentation.
  11 // Binaries may be compiled with this software without any royalties or
  12 // restrictions.
  13 //
  14 // The University of North Carolina at Chapel Hill makes no representations
  15 // about the suitability of this software for any purpose. It is provided
  16 // "as is" without express or implied warranty.
  17
  18 //============================================================================
  19 // tri.cpp : triangle defs and intersection routines.
  20 //============================================================================
  21
  22 #include <math.h>
  23
  24 //----------------------------------------------------------------------------
  25 // CHECKS IF 2D POINT P IS IN A 2D TRI ABC.
  26 //----------------------------------------------------------------------------
  27 int PtInTri2D(float Px, float Py,
  28               float Ax, float Ay, float Bx, float By, float Cx, float Cy)
  29 {
  30   float dABx=Bx-Ax, dABy=By-Ay, dBCx=Cx-Bx, dBCy=Cy-By; // "REPEATS"
  31   if (dABx*dBCy < dABy*dBCx) // CW
  32   {
  33     if (dABx*(Py-Ay) >= dABy*(Px-Ax)) return(0);        // ABxAP
  34     if (dBCx*(Py-By) >= dBCy*(Px-Bx)) return(0);        // BCxBP
  35     if ((Ax-Cx)*(Py-Cy) >= (Ay-Cy)*(Px-Cx)) return(0);  // CAxCP
  36   }
  37   else // CCW
  38   {
  39     if (dABx*(Py-Ay) < dABy*(Px-Ax)) return(0);         // ABxAP
  40     if (dBCx*(Py-By) < dBCy*(Px-Bx)) return(0);         // BCxBP
  41     if ((Ax-Cx)*(Py-Cy) < (Ay-Cy)*(Px-Cx)) return(0);   // CAxCP
  42   }
  43   return(1); // "INSIDE" EACH EDGE'S IN-HALF-SPACE (PT P IS INSIDE TRIANGLE)
  44 };
  45
  46 //----------------------------------------------------------------------------
  47 // CHECKS IF 3D POINT P (ON ABC'S PLANE) IS IN 3D TRI ABC.
  48 // P=PtOnTriPlane, N=PlaneNormal (does not have to be normalized)
  49 //----------------------------------------------------------------------------
  50 int PtInTri3D(const float P[3], float N[3],
  51               const float A[3], const float B[3], const float C[3])
  52 {
  53   #define ABS(x) (((x)<0)?(-(x)):x)  // HANDY UNIVERSAL ABSOLUTE VALUE FUNC
  54
  55   // DETERMINE LARGEST COMPONENT OF NORMAL (magnitude, since we want the largest projection)
  56   N[0]=ABS(N[0]); N[1]=ABS(N[1]); N[2]=ABS(N[2]);
  57
  58   // PROJECT ONTO PLANE WHERE PERPENDICULAR TO LARGEST NORMAL COMPONENT AXIS
  59   if (N[0]>N[1] && N[0]>N[2])      // X IS LARGEST SO PROJECT ONTO YZ-PLANE
  60     return( PtInTri2D(P[1],P[2], A[1],A[2], B[1],B[2], C[1],C[2]) );
  61   else if (N[1]>N[0] && N[1]>N[2]) // Y IS LARGEST SO PROJECT ONTO XZ-PLANE
  62     return( PtInTri2D(P[0],P[2], A[0],A[2], B[0],B[2], C[0],C[2]) );
  63   else                         // Z IS LARGEST SO PROJECT ONTO XY-PLANE
  64     return( PtInTri2D(P[0],P[1], A[0],A[1], B[0],B[1], C[0],C[1]) );
  65 }
  66
  67 //--------------------------------------------------------------------------
  68 // Checks if an edge UV intersects a triangle ABC. Returns whether or
  69 // not the edge intersects the plane (0 or 1) and the IsectPt.
  70 //--------------------------------------------------------------------------
  71 int EdgeTriIsect(const float U[3], const float V[3],
  72                  const float A[3], const float B[3], const float C[3],
  73                  float IsectPt[3])
  74 {
  75   // CALCULATE PLANE EQUATION COEFFICIENTS P FOR TRI ABC
  76   float P[4];
  77   float u[3] = {B[0]-A[0],B[1]-A[1],B[2]-A[2]};
  78   float v[3] = {C[0]-A[0],C[1]-A[1],C[2]-A[2]};
  79   P[0] = u[1]*v[2] - u[2]*v[1];  // CROSS-PROD BETWEEN u AND v
  80   P[1] = u[2]*v[0] - u[0]*v[2];  // DEFINES UNNORMALIZED NORMAL
  81   P[2] = u[0]*v[1] - u[1]*v[0];
  82   float l = (float)sqrt(P[0]*P[0] + P[1]*P[1] + P[2]*P[2]);  // NORMALIZE NORMAL
  83   P[0]/=l; P[1]/=l; P[2]/=l;
  84   P[3] = -(P[0]*A[0] + P[1]*A[1] + P[2]*A[2]);
  85
  86   // FIND INTERSECTION OF EDGE UV WITH PLANE P
  87   int EdgeIsectsPlane=0;
  88   float Dir[3] = { V[0]-U[0], V[1]-U[1], V[2]-U[2] };
  89   float NdotDir = P[0]*Dir[0] + P[1]*Dir[1] + P[2]*Dir[2];
  90   float NdotU = P[0]*U[0] + P[1]*U[1] + P[2]*U[2];
  91   if (NdotDir==(float)0) return(0);
  92   float t = (-P[3] - NdotU) / NdotDir;
  93   if (t>=0 && t<=1)
  94   {
  95     IsectPt[0] = U[0] + Dir[0]*t;
  96     IsectPt[1] = U[1] + Dir[1]*t;
  97     IsectPt[2] = U[2] + Dir[2]*t;
  98     EdgeIsectsPlane=1;
  99   }
 100
 101   // FIRST, DOES THE EDGE INTERSECT THE PLANE?
 102   if (!EdgeIsectsPlane) return(0);
 103
 104   // SEE IF ISECT PT IS IN THE TRI
 105   return( PtInTri3D(IsectPt,P,A,B,C) );
 106 }