home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / GRAPHICS / rayshade.lzh / poly.c < prev    next >
Text File  |  1990-05-08  |  8KB  |  310 lines
  1. /*
  2.  * poly.c
  3.  *
  4.  * Copyright (C) 1989, Craig E. Kolb
  5.  *
  6.  * This software may be freely copied, modified, and redistributed,
  7.  * provided that this copyright notice is preserved on all copies.
  8.  *
  9.  * There is no warranty or other guarantee of fitness for this software,
  10.  * it is provided solely .  Bug reports or fixes may be sent
  11.  * to the author, who may or may not act on them as he desires.
  12.  *
  13.  * You may not include this software in a program or other software product
  14.  * without supplying the source, or without informing the end-user that the
  15.  * source is available for no extra charge.
  16.  *
  17.  * If you modify this software, you should include a notice giving the
  18.  * name of the person performing the modification, the date of modification,
  19.  * and the reason for such modification.
  20.  *
  21.  * $Id: poly.c,v 3.0.1.2 90/04/04 19:01:19 craig Exp $
  22.  *
  23.  * $Log:    poly.c,v $
  24.  * Revision 3.0.1.2  90/04/04  19:01:19  craig
  25.  * patch5: Corrected free()-related problems.
  26.  * 
  27.  * Revision 3.0.1.1  89/12/06  16:33:38  craig
  28.  * patch2: Added calls to new error/warning routines.
  29.  * 
  30.  * Revision 3.0  89/10/27  02:05:59  craig
  31.  * Baseline for first official release.
  32.  * 
  33.  */
  34. #include <stdio.h>
  35. #include <math.h>
  36. #include "constants.h"
  37. #include "typedefs.h"
  38. #include "funcdefs.h"
  39.  
  40. /*
  41.  * Create a reference to a polygon with vertices equal to those
  42.  * on the linked-list "plist."
  43.  */
  44. Object *
  45. makpoly(surf, plist, npoints)
  46. char *surf;
  47. PointList *plist;
  48. int npoints;
  49. {
  50.     Polygon *poly;
  51.     Primitive *prim;
  52.     Object *newobj;
  53.     double indexval;
  54.     Vector edge1, edge2, anorm;
  55.     PointList *cur, *pltmp;
  56.     int i;
  57.     extern int TrashBadPoly;
  58.  
  59.     prim = mallocprim();
  60.     prim->type = POLY;
  61.     prim->surf = find_surface(surf);
  62.     poly = (Polygon *)Malloc(sizeof(Polygon));
  63.     prim->objpnt.p_poly = poly;
  64.     newobj = new_object(NULL, POLY, (char *)prim, (Trans *)NULL);
  65.     /*
  66.      * Allocate space for the vertices.
  67.      */
  68.     poly->points = (Vector *)Malloc((unsigned)(npoints*sizeof(Vector)));
  69.     poly->npoints = npoints;
  70.  
  71.     /*
  72.      * Copy the vertices from the linked list to the array, freeing
  73.      * the linked list as we go so that the caller doesn't have
  74.      * to worry about doing so.
  75.      */
  76.     i = npoints -1;
  77.     for(cur = plist; cur != (PointList *)0; cur = pltmp) {
  78.         poly->points[i--] = cur->vec;
  79.         pltmp = cur->next;
  80.         free((char *)cur);
  81.     }
  82.  
  83.     /*
  84.      * Find normal to polygon.  Check all edges before giving
  85.      * up, just to be relatively nice about things.
  86.      */
  87.     vecsub(poly->points[1], poly->points[0], &edge1);
  88.     for(i = 1;i < poly->npoints;i++) {
  89.         if(dotp(&edge1, &edge1) == 0.) {
  90.             if (TrashBadPoly) {
  91.                 free((char *)poly->points);
  92.                 free((char *)poly);
  93.                 free((char *)prim);
  94.                 free((char *)newobj);
  95.                 return (Object *)0;
  96.             }
  97.         }
  98.         vecsub(poly->points[(i+1)%npoints], poly->points[i], &edge2);
  99.         if(crossp(&poly->norm, &edge1, &edge2) != 0.)
  100.             break;
  101.         edge1 = edge2;
  102.     }
  103.  
  104.     if(i >= poly->npoints) {
  105.         /*
  106.           * If we walked all the way through the list,
  107.          * then we didn't find a valid normal vector -- we
  108.          * must have a degenerate polygon of some sort.
  109.          */
  110.         yywarning("Degenerate polygon.\n");
  111.         free((char *)poly->points);
  112.         free((char *)poly);
  113.         free((char *)prim);
  114.         free((char *)newobj);
  115.         return (Object *)0;
  116.     }
  117.  
  118.     /*
  119.      * Compute and store the plane constant.
  120.      */
  121.     poly->d = dotp(&poly->norm, &poly->points[0]);
  122.  
  123.     /*
  124.      * Find which part of the normal vector is "dominant."  This
  125.      * is used to turn the point-in-polygon test into a 2D problem.
  126.      */
  127.     anorm.x = abs(poly->norm.x);
  128.     anorm.y = abs(poly->norm.y);
  129.     anorm.z = abs(poly->norm.z);
  130.     indexval = max(anorm.y, anorm.z);
  131.     indexval = max(anorm.x, indexval);
  132.  
  133.     if(indexval == anorm.x)
  134.         poly->index = XNORMAL;
  135.     else if(indexval == anorm.y)
  136.         poly->index = YNORMAL;
  137.     else
  138.         poly->index = ZNORMAL;
  139.  
  140.     return newobj;
  141. }
  142.  
  143. /*
  144.  * Quadrants are defined as:
  145.  *        |
  146.  *   1    |   0
  147.  *        |
  148.  * -------c--------
  149.  *        |
  150.  *   2    |   3
  151.  *        |
  152.  */
  153. #define quadrant(p, c) ((p.u<c.u) ? ((p.v<c.v) ? 2 : 1) : ((p.v<c.v) ? 3 : 0))
  154.  
  155. /*
  156.  * Project a point in 3-space to the plane whose normal is indicated by "i."
  157.  */
  158. #define project(r, p, i)    {switch(i) { \
  159.                 case XNORMAL: \
  160.                     r.u = p.y; \
  161.                     r.v = p.z; \
  162.                     break; \
  163.                 case YNORMAL: \
  164.                     r.u = p.x; \
  165.                     r.v = p.z; \
  166.                     break; \
  167.                 case ZNORMAL: \
  168.                     r.u = p.x; \
  169.                     r.v = p.y; \
  170.                     break; \
  171.                   } }
  172. /*
  173.  * Perform ray-polygon intersection test.
  174.  */
  175. double
  176. intpoly(pos, ray, obj)
  177. Vector *pos, *ray;
  178. Primitive *obj;
  179. {
  180.     register Polygon *poly;
  181.     register int winding, i;
  182.     int quad, lastquad;
  183.     double dist, left, right;
  184.     Vec2d center, cur, last;
  185.     extern unsigned long primtests[];
  186.  
  187.     primtests[POLY]++;
  188.     poly = obj->objpnt.p_poly;
  189.     /*
  190.      * First, find where ray hits polygon plane, projecting
  191.      * along the polygon's dominant normal component.
  192.      */
  193.  
  194.     dist = dotp(&poly->norm, ray);
  195.     if(dist == 0.)
  196.         /*
  197.           * No intersection with polygon plane.
  198.          */
  199.         return 0.;
  200.  
  201.     dist = (poly->d - dotp(&poly->norm, pos)) / dist;
  202.     if(dist <= 0.)
  203.         /*
  204.          * The intersection point is behind the ray origin.
  205.          */
  206.         return 0.;
  207.  
  208.     /*
  209.      * Compute the point of intersection, projected appropriately.
  210.      */
  211.     if(poly->index == XNORMAL) {
  212.         center.u = pos->y + dist * ray->y;
  213.         center.v = pos->z + dist * ray->z;
  214.     } else if(poly->index == YNORMAL) {
  215.         center.v = pos->z + dist * ray->z;
  216.         center.u = pos->x + dist * ray->x;
  217.     } else {
  218.         center.u = pos->x + dist * ray->x;
  219.         center.v = pos->y + dist * ray->y;
  220.     }
  221.  
  222.     /*
  223.      * Is the point inside the polygon?
  224.      *
  225.      * Compute the winding number by finding the quadrant each
  226.      * polygon point lies in with respect to the the point in
  227.      * question, and computing a "delta" (winding number).  If we
  228.      * end up going around in a complete circle around
  229.      * the point (winding number is non-zero at the end), then
  230.      * we're inside.  Otherwise, the point is outside.
  231.      *
  232.      * Note that we can turn this into a 2D problem by projecting
  233.      * all the points along the axis defined by poly->index, which
  234.      * is the "dominant" part of the polygon's normal vector.
  235.      */
  236.     winding = 0;
  237.     project(last, poly->points[poly->npoints -1], poly->index);
  238.     lastquad = quadrant(last, center);
  239.     for(i = 0;i < poly->npoints; i++) {
  240.         project(cur, poly->points[i], poly->index);
  241.         quad = quadrant(cur, center);
  242.         if(quad != lastquad) {
  243.             if(((lastquad + 1) & 3) == quad)
  244.                 winding++;
  245.             else if(((quad + 1) & 3) == lastquad)
  246.                 winding--;
  247.             else {
  248.                 /*
  249.                  * Find where edge crosses
  250.                  * center's X axis.
  251.                  */
  252.                 right = last.u - cur.u;
  253.                 left = last.v - cur.v;
  254.                 left *= center.u - last.u;
  255.                 if(left + last.v * right > right * center.v)
  256.                     winding += 2;
  257.                 else
  258.                     winding -= 2;
  259.             }
  260.             lastquad = quad;
  261.         }
  262.         last = cur;
  263.  
  264.     }
  265.     return (winding == 0 ? 0. : dist);
  266. }
  267.  
  268. /*
  269.  * Return the normal to the polygon surface.
  270.  */
  271. /*ARGSUSED*/
  272. nrmpoly(pos, obj, nrm)
  273. Vector *pos, *nrm;
  274. Primitive *obj;
  275. {
  276.     *nrm = obj->objpnt.p_poly->norm;
  277. }
  278.  
  279. /*
  280.  * Compute the extent of a polygon
  281.  */
  282. polyextent(obj, bounds)
  283. Primitive *obj;
  284. double bounds[2][3];
  285. {
  286.     register Polygon *poly;
  287.     register int i;
  288.  
  289.     poly = obj->objpnt.p_poly;
  290.  
  291.     bounds[LOW][X] = bounds[HIGH][X] = poly->points[0].x;
  292.     bounds[LOW][Y] = bounds[HIGH][Y] = poly->points[0].y;
  293.     bounds[LOW][Z] = bounds[HIGH][Z] = poly->points[0].z;
  294.  
  295.     for(i = 1;i < poly->npoints;i++) {
  296.         if(poly->points[i].x < bounds[LOW][X])
  297.             bounds[LOW][X] = poly->points[i].x;
  298.         if(poly->points[i].x > bounds[HIGH][X])
  299.             bounds[HIGH][X] = poly->points[i].x;
  300.         if(poly->points[i].y < bounds[LOW][Y])
  301.             bounds[LOW][Y] = poly->points[i].y;
  302.         if(poly->points[i].y > bounds[HIGH][Y])
  303.             bounds[HIGH][Y] = poly->points[i].y;
  304.         if(poly->points[i].z < bounds[LOW][Z])
  305.             bounds[LOW][Z] = poly->points[i].z;
  306.         if(poly->points[i].z > bounds[HIGH][Z])
  307.             bounds[HIGH][Z] = poly->points[i].z;
  308.     }
  309. }
  310.