home *** CD-ROM | disk | FTP | other *** search
/ Compendium Deluxe 2 / LSD and 17bit Compendium Deluxe - Volume II.iso / a / prog / cprog / illum.lha / geo.c < prev    next >
Encoding:
C/C++ Source or Header  |  1988-10-21  |  7.3 KB  |  224 lines
  1. /* ****************************************************************
  2.  *                            geo.c
  3.  * ****************************************************************
  4.  *  MODULE PURPOSE:
  5.  *      This module contains routines that preform common
  6.  *      geometric manipulations for image synthesis.
  7.  *
  8.  *  MODULE CONTENTS:
  9.  *      geo_dot         - vector dot product
  10.  *      geo_cross       - vector cross product
  11.  *      geo_norm        - vector normalization
  12.  *      geo_line        - generate a line
  13.  *      geo_rfl         - compute the reflected vector
  14.  *      geo_rfr         - compute the refracted vector
  15.  *      geo_H           - compute H vector
  16.  *      geo_Ht          - compute the H' vector
  17.  */
  18. #include <stdio.h>
  19. #include <math.h>
  20. #include "geo.h"
  21.  
  22. /* ****************************************************************
  23.  *  double geo_dot (v0, v1)
  24.  *  DIR_VECT    *v0, *v1    (in) normalized direction vectors
  25.  *
  26.  *  Returns the dot product of two direction vectors, ref
  27.  *   Eq.(\*(rK).  Note that the dot product is the cosine
  28.  *   between the direction vectors because the direction
  29.  *   vectors are normalized.
  30.  */
  31. double geo_dot (v0, v1)
  32. DIR_VECT        *v0, *v1;
  33. {   return (v0->i * v1->i) + (v0->j * v1->j) + (v0->k * v1->k);
  34. }
  35.  
  36. /* ****************************************************************
  37.  *  DIR_VECT geo_cross (v0, v1)
  38.  *   DIR_VECT   *v0, *v1    (in)  normalized direction vectors
  39.  *
  40.  *  Returns the cross product of two direction vectors, refer
  41.  *   to Eq.(\*(rL).
  42.  */
  43. DIR_VECT geo_cross (v0, v1)
  44. DIR_VECT        *v0, *v1;
  45. {   DIR_VECT        v2;
  46.     v2.i = (v0->j * v1->k) - (v0->k * v1->j);
  47.     v2.j = (v0->k * v1->i) - (v0->i * v1->k);
  48.     v2.k = (v0->i * v1->j) - (v0->j * v1->i);
  49.     return v2;
  50. }
  51.  
  52. /* ****************************************************************
  53.  *  double geo_norm (*v0)
  54.  *   DIR_VECT   *v0         (mod) vector to be normalized
  55.  *
  56.  *  Returns the length of the vector before normalization.
  57.  *   Returns zero if the vector could not be normalized.
  58.  *   Note that the input vector is modified.
  59.  */
  60. double geo_norm (v0)
  61. DIR_VECT        *v0;
  62. {   double          len;
  63.     if ((len = geo_dot(v0, v0)) <= 0.0) return FALSE;
  64.     len = sqrt((double)len);
  65.     v0->i /= len;
  66.     v0->j /= len;
  67.     v0->k /= len;
  68.     return TRUE;
  69. }
  70.  
  71. /* ****************************************************************
  72.  *  double geo_line (*p0, *p1, *v)
  73.  *   POINT3     *p0, *p1    (in)  from, to points
  74.  *   LINE       *v          (out) generated line
  75.  *
  76.  *  Returns the distance from p0 to p1.  Returns 0.0
  77.  *   on error (p0=p1).  The line is expressed in parametric
  78.  *   form with a start point (p0) and a normalized direction
  79.  *   vector. The vector direction of the line is normalized.
  80.  */
  81. double geo_line (p0, p1, v)
  82. POINT3          *p0, *p1;
  83. LINE            *v;
  84. {   v->start = *p0;
  85.     v->dir.i = p1->x - p0->x;
  86.     v->dir.j = p1->y - p0->y;
  87.     v->dir.k = p1->z - p0->z;
  88.     return geo_norm(&v->dir);
  89. }
  90.  
  91. /* ****************************************************************
  92.  *  DIR_VECT *geo_rfl (V, N)
  93.  *   DIR_VECT   *V          (in)  incident vector
  94.  *   DIR_VECT   *N          (in)  surface normal
  95.  *
  96.  *  Returns the reflected direction vector.  The reflected
  97.  *   direction is computed using the method given by Whitted
  98.  *   (1980), refer to Eq.(\*(rM).
  99.  */
  100. DIR_VECT *geo_rfl (V, N)
  101. DIR_VECT        *V;
  102. DIR_VECT        *N;
  103. {   double              N_dot_V;
  104.     static DIR_VECT     rfl;
  105.     N_dot_V = geo_dot (N,V);
  106.     rfl.i = (2.0 * N_dot_V * N->i) - V->i;
  107.     rfl.j = (2.0 * N_dot_V * N->j) - V->j;
  108.     rfl.k = (2.0 * N_dot_V * N->k) - V->k;
  109.     return &rfl;
  110. }
  111.  
  112. /* ****************************************************************
  113.  *  DIR_VECT *geo_rfr (V, N, ni, nt)
  114.  *   DIR_VECT   *V          (in)  incident vector
  115.  *   DIR_VECT   *N          (in)  surface normal
  116.  *   double     ni          (in)  index of refraction for the
  117.  *                              material on the front of the
  118.  *                              interface (same side as N)
  119.  *   double     nt          (in)  index of refraction for the
  120.  *                              material on the back of the
  121.  *                              interface (opposite size as N)
  122.  *
  123.  *  Returns the refracted vector, if there complete internal
  124.  *   refracted (no refracted vector) then a NULL vector is
  125.  *   NULL is returned.  The vector is computed using the
  126.  *   method given by Hall (1983) with enhancements as
  127.  *   described in Appendix III.
  128.  */
  129. DIR_VECT *geo_rfr (V, N, ni, nt)
  130. DIR_VECT        *V;
  131. DIR_VECT        *N;
  132. double          ni;
  133. double          nt;
  134. {   static DIR_VECT T;      /* the refracted vector */
  135.     DIR_VECT    sin_T;      /* sin vect of the refracted vect */
  136.     DIR_VECT    cos_V;      /* cos vect of the incident vect */
  137.     double      len_sin_T;  /* length of sin T squared */
  138.     double      n_mult;     /* ni over nt */
  139.     double      N_dot_V;
  140.     double      N_dot_T;
  141.     if ((N_dot_V = geo_dot(N,V)) > 0.0) n_mult = ni / nt;
  142.     else n_mult = nt / ni;
  143.     cos_V.i = N_dot_V * N->i;
  144.     cos_V.j = N_dot_V * N->j;
  145.     cos_V.k = N_dot_V * N->k;
  146.     sin_T.i = (n_mult) * (cos_V.i - V->i);
  147.     sin_T.j = (n_mult) * (cos_V.j - V->j);
  148.     sin_T.k = (n_mult) * (cos_V.k - V->k);
  149.     if ((len_sin_T = geo_dot(&sin_T, &sin_T)) >= 1.0)
  150.     return NULL;    /* internal reflection */
  151.     N_dot_T = sqrt(1.0 - len_sin_T);
  152.     if (N_dot_V < 0.0) N_dot_T = -N_dot_T;
  153.     T.i = sin_T.i - (N->i * N_dot_T);
  154.     T.j = sin_T.j - (N->j * N_dot_T);
  155.     T.k = sin_T.k - (N->k * N_dot_T);
  156.     return &T;
  157. }
  158.  
  159. /* ****************************************************************
  160.  *  DIR_VECT *geo_H (L, V)
  161.  *   DIR_VECT   *L          (in) incident vector
  162.  *   DIR_VECT   *V          (in) reflected vector
  163.  *
  164.  *  Returns H, NULL on error (if L+H = 0).  Refer
  165.  *   to Eq.(\*(rO).
  166.  */
  167. DIR_VECT *geo_H (L, V)
  168. DIR_VECT        *L;
  169. DIR_VECT        *V;
  170. {   static DIR_VECT         H;
  171.     H.i = L->i + V->i;
  172.     H.j = L->j + V->j;
  173.     H.k = L->k + V->k;
  174.     if (!geo_norm(&H)) return NULL;
  175.     return &H;
  176. }
  177.  
  178. /* ****************************************************************
  179.  *  DIR_VECT *geo_Ht(L, T, ni, nt)
  180.  *   DIR_VECT   *L          (in) incident vector
  181.  *   DIR_VECT   *T          (in) transmitted hector
  182.  *   double     ni          (in) incident index
  183.  *   double     nt          (in) transmitted index
  184.  *
  185.  *  Returns H' oriented to the same side of the surface
  186.  *   as L computed using the method suggested by
  187.  *   Hall (1983).   Returns NULL on error (if the angle
  188.  *   between V and L is less than the critical angle).
  189.  */
  190. DIR_VECT *geo_Ht (L, T, ni, nt)
  191. DIR_VECT        *L;
  192. DIR_VECT        *T;
  193. double          ni;
  194. double          nt;
  195. {
  196.     float                   L_dot_T;
  197.     float                   divisor;
  198.     static DIR_VECT         Ht;
  199.  
  200.     L_dot_T = -(geo_dot(L,T));
  201.     /* check for special cases */
  202.     if (ni == nt) {
  203.     if (L_dot_T == 1.0) return L;
  204.     else return NULL;
  205.     }
  206.     if (ni < nt) {
  207.     if (L_dot_T < ni/nt) return NULL;
  208.     divisor = (nt / ni) - 1.0;
  209.     Ht.i = -(((L->i + T->i) / divisor) + T->i);
  210.     Ht.j = -(((L->j + T->j) / divisor) + T->j);
  211.     Ht.k = -(((L->k + T->k) / divisor) + T->k);
  212.     }
  213.     else
  214.     {   if (L_dot_T < nt/ni) return NULL;
  215.     divisor = (ni / nt) - 1.0;
  216.     Ht.i = ((L->i + T->i) / divisor) + L->i;
  217.     Ht.j = ((L->j + T->j) / divisor) + L->j;
  218.     Ht.k = ((L->k + T->k) / divisor) + L->k;
  219.     }
  220.     (void)geo_norm(&Ht);
  221.     return &Ht;
  222. }
  223. /* ************************************************************* */
  224.