--- ray/src/common/fvect.c 1991/03/19 13:14:19 1.5 +++ ray/src/common/fvect.c 2003/02/22 02:07:22 2.6 @@ -1,20 +1,71 @@ -/* Copyright (c) 1986 Regents of the University of California */ - #ifndef lint -static char SCCSid[] = "$SunId$ LBL"; +static const char RCSid[] = "$Id: fvect.c,v 2.6 2003/02/22 02:07:22 greg Exp $"; #endif - /* - * fvect.c - routines for float vector calculations + * fvect.c - routines for floating-point vector calculations + */ + +/* ==================================================================== + * The Radiance Software License, Version 1.0 * - * 8/14/85 + * Copyright (c) 1990 - 2002 The Regents of the University of California, + * through Lawrence Berkeley National Laboratory. All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in + * the documentation and/or other materials provided with the + * distribution. + * + * 3. The end-user documentation included with the redistribution, + * if any, must include the following acknowledgment: + * "This product includes Radiance software + * (http://radsite.lbl.gov/) + * developed by the Lawrence Berkeley National Laboratory + * (http://www.lbl.gov/)." + * Alternately, this acknowledgment may appear in the software itself, + * if and wherever such third-party acknowledgments normally appear. + * + * 4. The names "Radiance," "Lawrence Berkeley National Laboratory" + * and "The Regents of the University of California" must + * not be used to endorse or promote products derived from this + * software without prior written permission. For written + * permission, please contact radiance@radsite.lbl.gov. + * + * 5. Products derived from this software may not be called "Radiance", + * nor may "Radiance" appear in their name, without prior written + * permission of Lawrence Berkeley National Laboratory. + * + * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED + * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR + * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF + * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT + * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * ==================================================================== + * + * This software consists of voluntary contributions made by many + * individuals on behalf of Lawrence Berkeley National Laboratory. For more + * information on Lawrence Berkeley National Laboratory, please see + * . */ +#include #include "fvect.h" -#define FTINY 1e-7 - double fdot(v1, v2) /* return the dot product of two vectors */ register FVECT v1, v2; @@ -27,11 +78,12 @@ double dist2(p1, p2) /* return square of distance between points */ register FVECT p1, p2; { - static FVECT delta; + FVECT delta; delta[0] = p2[0] - p1[0]; delta[1] = p2[1] - p1[1]; delta[2] = p2[2] - p1[2]; + return(DOT(delta, delta)); } @@ -41,13 +93,13 @@ dist2line(p, ep1, ep2) /* return square of distance t FVECT p; /* the point */ FVECT ep1, ep2; /* points on the line */ { - static double d, d1, d2; + register double d, d1, d2; d = dist2(ep1, ep2); d1 = dist2(ep1, p); - d2 = dist2(ep2, p); + d2 = d + d1 - dist2(ep2, p); - return(d1 - (d+d1-d2)*(d+d1-d2)/d/4); + return(d1 - 0.25*d2*d2/d); } @@ -56,7 +108,7 @@ dist2lseg(p, ep1, ep2) /* return square of distance t FVECT p; /* the point */ FVECT ep1, ep2; /* the end points */ { - static double d, d1, d2; + register double d, d1, d2; d = dist2(ep1, ep2); d1 = dist2(ep1, p); @@ -69,11 +121,13 @@ FVECT ep1, ep2; /* the end points */ if (d1 - d2 > d) return(d2); } + d2 = d + d1 - d2; - return(d1 - (d+d1-d2)*(d+d1-d2)/d/4); /* distance to line */ + return(d1 - 0.25*d2*d2/d); /* distance to line */ } +void fcross(vres, v1, v2) /* vres = v1 X v2 */ register FVECT vres, v1, v2; { @@ -83,9 +137,10 @@ register FVECT vres, v1, v2; } +void fvsum(vres, v0, v1, f) /* vres = v0 + f*v1 */ -FVECT vres, v0, v1; -double f; +register FVECT vres, v0, v1; +register double f; { vres[0] = v0[0] + f*v1[0]; vres[1] = v0[1] + f*v1[1]; @@ -97,45 +152,44 @@ double normalize(v) /* normalize a vector, return old magnitude */ register FVECT v; { - static double len; + register double len, d; - len = DOT(v, v); + d = DOT(v, v); - if (len <= 0.0) + if (d <= 0.0) return(0.0); - /****** problematic - if (len >= (1.0-FTINY)*(1.0-FTINY) && - len <= (1.0+FTINY)*(1.0+FTINY)) - return(1.0); - ******/ + if (d <= 1.0+FTINY && d >= 1.0-FTINY) + len = 0.5 + 0.5*d; /* first order approximation */ + else + len = sqrt(d); - len = sqrt(len); - v[0] /= len; - v[1] /= len; - v[2] /= len; + v[0] *= d = 1.0/len; + v[1] *= d; + v[2] *= d; + return(len); } +void spinvector(vres, vorig, vnorm, theta) /* rotate vector around normal */ FVECT vres, vorig, vnorm; double theta; { - extern double sin(), cos(); - double sint, cost, dotp; + double sint, cost, normprod; FVECT vperp; register int i; if (theta == 0.0) { - VCOPY(vres, vorig); + if (vres != vorig) + VCOPY(vres, vorig); return; } - sint = sin(theta); cost = cos(theta); - dotp = DOT(vorig, vnorm); + sint = sin(theta); + normprod = DOT(vorig, vnorm)*(1.-cost); fcross(vperp, vnorm, vorig); for (i = 0; i < 3; i++) - vres[i] = vnorm[i]*dotp*(1.-cost) + - vorig[i]*cost + vperp[i]*sint; + vres[i] = vorig[i]*cost + vnorm[i]*normprod + vperp[i]*sint; }