VecUtil.cc

00001 /*
00002 File: VecUtil.cc
00003 
00004 Function: Contains a grab-bag of useful graphics
00005 vector routines.
00006 
00007 Author(s): Andrew Willmott
00008 
00009 */
00010 
00011 #include "gcl/VecUtil.h"
00012 
00013  Bool Refract(
00014 GCLReal fromIndex, 
00015 GCLReal toIndex,
00016 const Vector &v, 
00017 const Vector &n, 
00018 GCLReal cosNV, 
00019 Vector &refractDir
00020 )
00022 {
00023 GCLReal kn, cos2, k;
00024 
00025 kn = fromIndex / toIndex;
00026 cos2 = 1.0 - sqr(kn) * (1.0 - sqr(cosNV));
00027 
00028 if (cos2 < 0.0)
00029 return(1); // Total internal reflection!
00030 
00031 k = kn * cosNV - sqrt(cos2);
00032 
00033 refractDir = -kn * v + k * n;
00034 
00035  return(0);
00036 }
00037 
00038 Void CalcTriAreaNormal(
00039 const Point &a,
00040 const Point &b, 
00041 const Point &c,
00042 Vector &n
00043 )
00048 {
00049 n[0] = (a[1] - b[1]) * (a[2] + b[2]) + (b[1] - c[1]) * (b[2] + c[2])
00050 + (c[1] - a[1]) * (c[2] + a[2]);
00051 n[1] = (a[2] - b[2]) * (a[0] + b[0]) + (b[2] - c[2]) * (b[0] + c[0])
00052 + (c[2] - a[2]) * (c[0] + a[0]);
00053 n[2] = (a[0] - b[0]) * (a[1] + b[1]) + (b[0] - c[0]) * (b[1] + c[1]) 
00054  + (c[0] - a[0]) * (c[1] + a[1]);
00055 }
00056 
00057 Void UpdateBounds(const Point& pt, Point& min, Point& max)
00058 {
00059 if (min[0] > pt[0])
00060 min[0] = pt[0];
00061 else if (max[0] < pt[0])
00062 max[0] = pt[0];
00063 
00064 if (min[1] > pt[1])
00065 min[1] = pt[1];
00066 else if (max[1] < pt[1])
00067 max[1] = pt[1];
00068 
00069 if (min[2] > pt[2])
00070 min[2] = pt[2];
00071 else if (max[2] < pt[2])
00072  max[2] = pt[2];
00073 }
00074 
00075  inline Double xCross(const Vector &a, const Vector &b, const Vector &c)
00076 // returns (a - b) x (c - a) . (1, 0, 0)
00077 { return((a[1] - b[1]) * (c[2] - b[2]) - (a[2] - b[2]) * (c[1] - b[1])); }
00078  inline Double yCross(const Vector &a, const Vector &b, const Vector &c)
00079 // returns (a - b) x (c - a) . (0, 1, 0)
00080 { return((a[2] - b[2]) * (c[0] - b[0]) - (a[0] - b[0]) * (c[2] - b[2])); }
00081 inline Double zCross(const Vector &a, const Vector &b, const Vector &c)
00082  // returns (a - b) x (c - a) . (0, 0, 1)
00083 { return((a[0] - b[0]) * (c[1] - b[1]) - (a[1] - b[1]) * (c[0] - b[0])); }
00084  
00085 const GCLReal kEdgeFuzz = 0.0;
00086 
00087 Bool PointIsInsideTriangle(
00088 const Point &p0, 
00089 const Point &p1, 
00090 const Point &p2, 
00091 const Point &point,
00092 Vector *coords 
00093 )
00097 {
00098 GCLReal temp, eps, a, b, c;
00099 Vector normal;
00100 Int normMajorAxis;
00101 
00102 // Calculate face normal & other info
00103 
00104 CalcTriAreaNormal(p0, p1, p2, normal);
00105 normal[0] = abs(normal[0]);
00106 normal[1] = abs(normal[1]);
00107 normal[2] = abs(normal[2]);
00108 normMajorAxis = MaxEltIndex(normal); 
00109 temp = normal[normMajorAxis];
00110 if (temp == 0.0)
00111 // degenerate triangle... can't hit this.
00112 return(false);
00113 
00114 eps = kEdgeFuzz * temp;
00115 switch(normMajorAxis)
00116 {
00117 case 0:
00118 a = xCross(p2, p1, point);
00119 if (a < eps) return(false);
00120 
00121 b = xCross(p0, p2, point);
00122 if (b < eps) return(false);
00123 
00124 c = temp - a - b;
00125 if (c < eps) return(false);
00126 
00127 break;
00128 
00129 case 1:
00130 a = yCross(p2, p1, point);
00131 if (a < eps) return(false);
00132 
00133 b = yCross(p0, p2, point); 
00134 if (b < eps) return(false);
00135 
00136 c = temp - a - b;
00137 if (c < eps) return(false);
00138 
00139 break;
00140 
00141 case 2:
00142 a = zCross(p2, p1, point);
00143 if (a < eps) return(false);
00144 
00145 b = zCross(p0, p2, point);
00146 if (b < eps) return(false);
00147 
00148 c = temp - a - b;
00149 if (c < eps) return(false);
00150 
00151 break;
00152 }
00153 
00154 if (coords)
00155 {
00156 (*coords)[0] = a / temp;
00157 (*coords)[1] = b / temp;
00158 (*coords)[2] = c / temp;
00159 }
00160 
00161  return(true);
00162 }
00163 
00164 VecTrans ZViewMatrix3(Vector &dir)
00167 {
00168 VecTrans M;
00169 
00170 M[0] = norm(MVector(dir[0], dir[1], dir[2]));
00171 M[1] = norm(MVector(dir[1], -dir[0], 0.0));
00172 M[2] = norm(MVector(-dir[0] * dir[2], -dir[1] * dir[2], 
00173 sqr(dir[0]) + sqr(dir[1])));
00174 
00175 #ifndef VL_ROW_ORIENT
00176 M = trans(M);
00177 #endif
00178 
00179  return(M);
00180 }
00181 
00182 Void UnitSquareToUnitDisc(Coord in, Coord &out)
00207 {
00208 GCLReal phi, r;
00209 Coord a = 2.0 * in - vl_1; /* a is now on [-1,1]^2 */
00210 
00211 if (a[0] > -a[1])
00212 /* region 1 or 2 */
00213 {
00214 if (a[0] > a[1])
00215 /* region 1, also |a| > |b| */
00216 {
00217 r = a[0];
00218 phi = (vl_pi / 4.0) * (a[1] / a[0]);
00219 }
00220 else
00221 /* region 2, also |b| > |a| */
00222 {
00223 r = a[1];
00224 phi = (vl_pi / 4.0) * (2.0 - (a[0] / a[1]));
00225 }
00226 }
00227 else 
00228 /* region 3 or 4 */
00229 {
00230 if (a[0] < a[1]) 
00231 /* region 3, also |a| >= |b|, a != 0 */
00232 { 
00233 r = -a[0];
00234 phi = (vl_pi / 4.0) * (4.0 + (a[1] / a[0]));
00235 }
00236 else 
00237 /* region 4, |b| >= |a|, but a==0 and b==0 could occur. */
00238 {
00239 r = -a[1];
00240 if (a[1] != 0)
00241 phi = (vl_pi / 4.0) * (6.0 - (a[0] / a[1]));
00242 else
00243 phi = 0.0;
00244 }
00245 }
00246 
00247 out[0] = r * cos(phi);
00248  out[1] = r * sin(phi);
00249 }
00250 
00251 Void UnitDiskToUnitSquare(Coord in, Coord &out)
00257 {
00258 GCLReal r = len(in);
00259 GCLReal phi = atan2(in[0], in[1]);
00260 Coord a;
00261 
00262 if (phi < -vl_pi / 4.0) 
00263 /* in range [-pi/4,7pi/4] */
00264 phi += 2.0 * vl_pi; 
00265 if (phi < vl_pi / 4.0)
00266 /* region 1 */
00267 { 
00268 a[0] = r;
00269 a[1] = phi * a[0] / (vl_pi / 4.0);
00270 }
00271 else if (phi < 3.0 * vl_pi / 4.0)
00272 /* region 2 */
00273 {
00274 a[1] = r;
00275 a[0] = -(phi - vl_pi / 2.0) * a[1] / (vl_pi / 4.0);
00276 }
00277 else if (phi < 5.0 * vl_pi / 4.0)
00278 /* region 3 */
00279 {
00280 a[0] = -r;
00281 a[1] = (phi - vl_pi) * a[0] / (vl_pi / 4.0);
00282 }
00283 else
00284 /* region 4 */
00285 {
00286 a[1] = -r;
00287 a[0] = -(phi - 3.0 * vl_pi / 2.0) * a[1] / (vl_pi / 4.0);
00288 }
00289 
00290 out[0] = (a[0] + 1) / 2;
00291  out[1] = (a[1] + 1) / 2;
00292 }
00293 
00294 GCLReal ProjectToBox(Point &min, Point &max, Point &p, Vector &r)
00297 {
00298 Vector d;
00299 Int i;
00300 
00301 for (i = 0; i < 3; i++)
00302 if (r[i] > 0.0)
00303 d[i] = (max[i] - p[i]) / r[i];
00304 else if (r[i] < 0.0)
00305 d[i] = (min[i] - p[i]) / r[i];
00306 else
00307 d[i] = HUGE_VAL;
00308 
00309  return(MinElt(d));
00310 }
00311 
00312 Vector RandomVector()
00313 {
00314 Vector v;
00315 
00316 v[0] = vl_rand();
00317 v[1] = vl_rand();
00318 v[2] = vl_rand();
00319 
00320 return(v);
00321 }