src/rt/srcdraw.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#endif
/*
 * Draw small sources into image in case we missed them.
 *
 *  External symbols declared in ray.h
 */

/* ====================================================================
 * The Radiance Software License, Version 1.0
 *
 * 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 [email protected].
 *
 * 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
 * <http://www.lbl.gov/>.
 */

#include  "ray.h"

#include  "view.h"

#include  "source.h"


#define CLIP_ABOVE      1
#define CLIP_BELOW      2
#define CLIP_RIGHT      3
#define CLIP_LEFT       4

#define MAXVERT         10

typedef struct splist {
        struct splist   *next;                  /* next source in list */
        int     sn;                             /* source number */
        short   nv;                             /* number of vertices */
        FLOAT   vl[3][2];                       /* vertex array (last) */
} SPLIST;                               /* source polygon list */

extern VIEW     ourview;                /* our view parameters */
extern int      hres, vres;             /* our image resolution */
static SPLIST   *sphead = NULL;         /* our list of source polys */


static int
inregion(p, cv, crit)                   /* check if vertex is in region */
FLOAT   p[2];
double  cv;
int     crit;
{
        switch (crit) {
        case CLIP_ABOVE:
                return(p[1] < cv);
        case CLIP_BELOW:
                return(p[1] >= cv);
        case CLIP_RIGHT:
                return(p[0] < cv);
        case CLIP_LEFT:
                return(p[0] >= cv);
        }
        return(-1);
}


static
clipregion(a, b, cv, crit, r)           /* find intersection with boundary */
register FLOAT  a[2], b[2];
double  cv;
int     crit;
FLOAT   r[2];   /* return value */
{
        switch (crit) {
        case CLIP_ABOVE:
        case CLIP_BELOW:
                r[1] = cv;
                r[0] = a[0] + (cv-a[1])/(b[1]-a[1])*(b[0]-a[0]);
                return;
        case CLIP_RIGHT:
        case CLIP_LEFT:
                r[0] = cv;
                r[1] = a[1] + (cv-a[0])/(b[0]-a[0])*(b[1]-a[1]);
                return;
        }
}


static int
hp_clip_poly(vl, nv, cv, crit, vlo)     /* clip polygon to half-plane */
FLOAT   vl[][2];
int     nv;
double  cv;
int     crit;
FLOAT   vlo[][2];       /* return value */
{
        FLOAT   *s, *p;
        register int    j, nvo;

        s = vl[nv-1];
        nvo = 0;
        for (j = 0; j < nv; j++) {
                p = vl[j];
                if (inregion(p, cv, crit)) {
                        if (!inregion(s, cv, crit))
                                clipregion(s, p, cv, crit, vlo[nvo++]);
                        vlo[nvo][0] = p[0]; vlo[nvo++][1] = p[1];
                } else if (inregion(s, cv, crit))
                        clipregion(s, p, cv, crit, vlo[nvo++]);
                s = p;
        }
        return(nvo);
}


static int
box_clip_poly(vl, nv, xl, xr, yb, ya, vlo)      /* clip polygon to box */
FLOAT   vl[MAXVERT][2];
int     nv;
double  xl, xr, yb, ya;
FLOAT   vlo[MAXVERT][2];        /* return value */
{
        FLOAT   vlt[MAXVERT][2];
        int     nvt, nvo;

        nvt = hp_clip_poly(vl, nv, yb, CLIP_BELOW, vlt);
        nvo = hp_clip_poly(vlt, nvt, ya, CLIP_ABOVE, vlo);
        nvt = hp_clip_poly(vlo, nvo, xl, CLIP_LEFT, vlt);
        nvo = hp_clip_poly(vlt, nvt, xr, CLIP_RIGHT, vlo);

        return(nvo);
}


static double
minw2(vl, nv, ar2)                      /* compute square of minimum width */
FLOAT   vl[][2];
int     nv;
double  ar2;
{
        double  d2, w2, w2min, w2max;
        register FLOAT  *p0, *p1, *p2;
        int     i, j;
                                /* find minimum for all widths */
        w2min = FHUGE;
        p0 = vl[nv-1];
        for (i = 0; i < nv; i++) {      /* for each edge */
                p1 = vl[i];
                d2 = (p1[0]-p0[0])*(p1[0]-p0[0]) +
                                (p1[1]-p0[1])*(p1[1]-p0[1])*ar2;
                w2max = 0.;             /* find maximum for this side */
                for (j = 1; j < nv-1; j++) {
                        p2 = vl[(i+j)%nv];
                        w2 = (p1[0]-p0[0])*(p2[1]-p0[1]) -
                                        (p1[1]-p0[1])*(p2[0]-p0[0]);
                        w2 = w2*w2*ar2/d2;      /* triangle height squared */
                        if (w2 > w2max)
                                w2max = w2;
                }
                if (w2max < w2min)      /* global min. based on local max.'s */
                        w2min = w2max;
                p0 = p1;
        }
        return(w2min);
}


static
convex_center(vl, nv, cv)               /* compute center of convex polygon */
register FLOAT  vl[][2];
int     nv;
FLOAT   cv[2];          /* return value */
{
        register int    i;
                                        /* simple average (suboptimal) */
        cv[0] = cv[1] = 0.;
        for (i = 0; i < nv; i++) {
                cv[0] += vl[i][0];
                cv[1] += vl[i][1];
        }
        cv[0] /= (double)nv;
        cv[1] /= (double)nv;
}


static double
poly_area(vl, nv)                       /* compute area of polygon */
register FLOAT  vl[][2];
int     nv;
{
        double  a;
        FLOAT   v0[2], v1[2];
        register int    i;

        a = 0.;
        v0[0] = vl[1][0] - vl[0][0];
        v0[1] = vl[1][1] - vl[0][1];
        for (i = 2; i < nv; i++) {
                v1[0] = vl[i][0] - vl[0][0];
                v1[1] = vl[i][1] - vl[0][1];
                a += v0[0]*v1[1] - v0[1]*v1[0];
                v0[0] = v1[0]; v0[1] = v1[1];
        }
        return(a * (a >= 0. ? .5 : -.5));
}


static int
convex_hull(vl, nv, vlo)                /* compute polygon's convex hull */
FLOAT   vl[][2];
int     nv;
FLOAT   vlo[][2];       /* return value */
{
        int     nvo, nvt;
        FLOAT   vlt[MAXVERT][2];
        double  voa, vta;
        register int    i, j;
                                        /* start with original polygon */
        for (i = nvo = nv; i--; ) {
                vlo[i][0] = vl[i][0]; vlo[i][1] = vl[i][1];
        }
        voa = poly_area(vlo, nvo);      /* compute its area */
        for (i = 0; i < nvo; i++) {             /* for each output vertex */
                for (j = 0; j < i; j++) {
                        vlt[j][0] = vlo[j][0]; vlt[j][1] = vlo[j][1];
                }
                nvt = nvo - 1;                  /* make poly w/o vertex */
                for (j = i; j < nvt; j++) {
                        vlt[j][0] = vlo[j+1][0]; vlt[j][1] = vlo[j+1][1];
                }
                vta = poly_area(vlt, nvt);
                if (vta >= voa) {               /* is simpler poly bigger? */
                        voa = vta;                      /* then use it */
                        for (j = nvo = nvt; j--; ) {
                                vlo[j][0] = vlt[j][0]; vlo[j][1] = vlt[j][1];
                        }
                        i--;                            /* next adjust */
                }
        }
        return(nvo);
}


static
spinsert(sn, vl, nv)                    /* insert new source polygon */
int     sn;
FLOAT   vl[][2];
int     nv;
{
        register SPLIST *spn;
        register int    i;

        if (nv < 3)
                return;
        if (nv > 3)
                spn = (SPLIST *)malloc(sizeof(SPLIST)+sizeof(FLOAT)*2*(nv-3));
        else
                spn = (SPLIST *)malloc(sizeof(SPLIST));
        if (spn == NULL)
                error(SYSTEM, "out of memory in spinsert");
        spn->sn = sn;
        for (i = spn->nv = nv; i--; ) {
                spn->vl[i][0] = vl[i][0]; spn->vl[i][1] = vl[i][1];
        }
        spn->next = sphead;             /* push onto global list */
        sphead = spn;
}


int
sourcepoly(sn, sp)                      /* compute image polygon for source */
int     sn;
FLOAT   sp[MAXVERT][2];
{
        static char     cubeord[8][6] = {{1,3,2,6,4,5},{0,4,5,7,3,2},
                                         {0,1,3,7,6,4},{0,1,5,7,6,2},
                                         {0,2,6,7,5,1},{0,4,6,7,3,1},
                                         {0,2,3,7,5,4},{1,5,4,6,2,3}};
        register SRCREC *s = source + sn;
        FVECT   ap, ip;
        FLOAT   pt[6][2];
        int     dir;
        register int    i, j;

        if (s->sflags & (SDISTANT|SFLAT)) {
                if (s->sflags & SDISTANT && ourview.type == VT_PAR)
                        return(0);              /* all or nothing case */
                if (s->sflags & SFLAT) {
                        for (i = 0; i < 3; i++)
                                ap[i] = s->sloc[i] - ourview.vp[i];
                        if (DOT(ap, s->snorm) >= 0.)
                                return(0);      /* source faces away */
                }
                for (j = 0; j < 4; j++) {       /* four corners */
                        for (i = 0; i < 3; i++) {
                                ap[i] = s->sloc[i];
                                if (j==1|j==2) ap[i] += s->ss[SU][i];
                                else ap[i] -= s->ss[SU][i];
                                if (j==2|j==3) ap[i] += s->ss[SV][i];
                                else ap[i] -= s->ss[SV][i];
                                if (s->sflags & SDISTANT) {
                                        ap[i] *= 1. + ourview.vfore;
                                        ap[i] += ourview.vp[i];
                                }
                        }
                        viewloc(ip, &ourview, ap);      /* find image point */
                        if (ip[2] <= 0.)
                                return(0);              /* in front of view */
                        sp[j][0] = ip[0]; sp[j][1] = ip[1];
                }
                return(4);
        }
                                        /* identify furthest corner */
        for (i = 0; i < 3; i++)
                ap[i] = s->sloc[i] - ourview.vp[i];
        dir =   (DOT(ap,s->ss[SU])>0.) |
                (DOT(ap,s->ss[SV])>0.)<<1 |
                (DOT(ap,s->ss[SW])>0.)<<2 ;
                                        /* order vertices based on this */
        for (j = 0; j < 6; j++) {
                for (i = 0; i < 3; i++) {
                        ap[i] = s->sloc[i];
                        if (cubeord[dir][j] & 1) ap[i] += s->ss[SU][i];
                        else ap[i] -= s->ss[SU][i];
                        if (cubeord[dir][j] & 2) ap[i] += s->ss[SV][i];
                        else ap[i] -= s->ss[SV][i];
                        if (cubeord[dir][j] & 4) ap[i] += s->ss[SW][i];
                        else ap[i] -= s->ss[SW][i];
                }
                viewloc(ip, &ourview, ap);      /* find image point */
                if (ip[2] <= 0.)
                        return(0);              /* in front of view */
                pt[j][0] = ip[0]; pt[j][1] = ip[1];
        }
        return(convex_hull(pt, 6, sp));         /* make sure it's convex */
}


                        /* initialize by finding sources smaller than rad */
init_drawsources(rad)
int     rad;                            /* source sample size */
{
        FLOAT   spoly[MAXVERT][2];
        int     nsv;
        register SPLIST *sp;
        register int    i;
                                        /* free old source list if one */
        for (sp = sphead; sp != NULL; sp = sphead) {
                sphead = sp->next;
                free((void *)sp);
        }
                                        /* loop through all sources */
        for (i = nsources; i--; ) {
                                        /* compute image polygon for source */
                if (!(nsv = sourcepoly(i, spoly)))
                        continue;
                                        /* clip to image boundaries */
                if (!(nsv = box_clip_poly(spoly, nsv, 0., 1., 0., 1., spoly)))
                        continue;
                                        /* big enough for standard sampling? */
                if (minw2(spoly, nsv, ourview.vn2/ourview.hn2) >
                                (double)rad*rad/hres/hres)
                        continue;
                                        /* OK, add to our list */
                spinsert(i, spoly, nsv);
        }
}

void                    /* add sources smaller than rad to computed subimage */
drawsources(pic, zbf, x0, xsiz, y0, ysiz)
COLOR   *pic[];                         /* subimage pixel value array */
float   *zbf[];                         /* subimage distance array (opt.) */
int     x0, xsiz, y0, ysiz;             /* origin and size of subimage */
{
        FLOAT   spoly[MAXVERT][2], ppoly[MAXVERT][2];
        int     nsv, npv;
        int     xmin, xmax, ymin, ymax, x, y;
        FLOAT   cxy[2];
        double  w;
        RAY     sr;
        register SPLIST *sp;
        register int    i;
                                        /* check each source in our list */
        for (sp = sphead; sp != NULL; sp = sp->next) {
                                        /* clip source poly to subimage */
                nsv = box_clip_poly(sp->vl, sp->nv,
                                (double)x0/hres, (double)(x0+xsiz)/hres,
                                (double)y0/vres, (double)(y0+ysiz)/vres, spoly);
                if (!nsv)
                        continue;
                                        /* find common subimage (BBox) */
                xmin = x0 + xsiz; xmax = x0;
                ymin = y0 + ysiz; ymax = y0;
                for (i = 0; i < nsv; i++) {
                        if ((double)xmin/hres > spoly[i][0])
                                xmin = spoly[i][0]*hres + FTINY;
                        if ((double)xmax/hres < spoly[i][0])
                                xmax = spoly[i][0]*hres - FTINY;
                        if ((double)ymin/vres > spoly[i][1])
                                ymin = spoly[i][1]*vres + FTINY;
                        if ((double)ymax/vres < spoly[i][1])
                                ymax = spoly[i][1]*vres - FTINY;
                }
                                        /* evaluate each pixel in BBox */
                for (y = ymin; y <= ymax; y++)
                        for (x = xmin; x <= xmax; x++) {
                                                        /* subarea for pixel */
                                npv = box_clip_poly(spoly, nsv,
                                                (double)x/hres, (x+1.)/hres,
                                                (double)y/vres, (y+1.)/vres,
                                                ppoly);
                                if (!npv)
                                        continue;       /* no overlap */
                                convex_center(ppoly, npv, cxy);
                                if ((sr.rmax = viewray(sr.rorg,sr.rdir,&ourview,
                                                cxy[0],cxy[1])) < -FTINY)
                                        continue;       /* not in view */
                                if (source[sp->sn].sflags & SSPOT &&
                                                spotout(&sr, source[sp->sn].sl.s))
                                        continue;       /* outside spot */
                                rayorigin(&sr, NULL, SHADOW, 1.0);
                                sr.rsrc = sp->sn;
                                rayvalue(&sr);          /* compute value */
                                if (bright(sr.rcol) <= FTINY)
                                        continue;       /* missed/blocked */
                                                        /* modify pixel */
                                if (zbf[y-y0] != NULL &&
                                                sr.rt < 0.999*zbf[y-y0][x-x0])
                                        zbf[y-y0][x-x0] = sr.rt;
                                else if (!bigdiff(sr.rcol, pic[y-y0][x-x0],
                                                0.001)) /* source sample */
                                        setcolor(pic[y-y0][x-x0], 0., 0., 0.);
                                w = poly_area(ppoly, npv) * hres * vres;
                                scalecolor(sr.rcol, w);
                                scalecolor(pic[y-y0][x-x0], 1.-w);
                                addcolor(pic[y-y0][x-x0], sr.rcol);
                        }
        }
}
Revision:	2.5
Committed:	Sat Feb 22 02:07:29 2003 UTC (21 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.4:	+62 -6 lines
Log Message:	Changes and check-in for 3.5 release Includes new source files and modifications not recorded for many years See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.5	static const char RCSid[] = "$Id$";
3	greg	2.1	#endif
4			/*
5			* Draw small sources into image in case we missed them.
6	greg	2.5	*
7			* External symbols declared in ray.h
8			*/
9
10			/* ====================================================================
11			* The Radiance Software License, Version 1.0
12			*
13			* Copyright (c) 1990 - 2002 The Regents of the University of California,
14			* through Lawrence Berkeley National Laboratory. All rights reserved.
15			*
16			* Redistribution and use in source and binary forms, with or without
17			* modification, are permitted provided that the following conditions
18			* are met:
19			*
20			* 1. Redistributions of source code must retain the above copyright
21			* notice, this list of conditions and the following disclaimer.
22			*
23			* 2. Redistributions in binary form must reproduce the above copyright
24			* notice, this list of conditions and the following disclaimer in
25			* the documentation and/or other materials provided with the
26			* distribution.
27			*
28			* 3. The end-user documentation included with the redistribution,
29			* if any, must include the following acknowledgment:
30			* "This product includes Radiance software
31			* (http://radsite.lbl.gov/)
32			* developed by the Lawrence Berkeley National Laboratory
33			* (http://www.lbl.gov/)."
34			* Alternately, this acknowledgment may appear in the software itself,
35			* if and wherever such third-party acknowledgments normally appear.
36			*
37			* 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
38			* and "The Regents of the University of California" must
39			* not be used to endorse or promote products derived from this
40			* software without prior written permission. For written
41			* permission, please contact [email protected].
42			*
43			* 5. Products derived from this software may not be called "Radiance",
44			* nor may "Radiance" appear in their name, without prior written
45			* permission of Lawrence Berkeley National Laboratory.
46			*
47			* THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
48			* WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
49			* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
50			* DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
51			* ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
52			* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
53			* LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
54			* USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
55			* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
56			* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
57			* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
58			* SUCH DAMAGE.
59			* ====================================================================
60			*
61			* This software consists of voluntary contributions made by many
62			* individuals on behalf of Lawrence Berkeley National Laboratory. For more
63			* information on Lawrence Berkeley National Laboratory, please see
64			* <http://www.lbl.gov/>.
65	greg	2.1	*/
66
67			#include "ray.h"
68
69			#include "view.h"
70
71			#include "source.h"
72
73
74			#define CLIP_ABOVE 1
75			#define CLIP_BELOW 2
76			#define CLIP_RIGHT 3
77			#define CLIP_LEFT 4
78
79			#define MAXVERT 10
80
81	greg	2.3	typedef struct splist {
82			struct splist next; / next source in list */
83			int sn; /* source number */
84			short nv; /* number of vertices */
85			FLOAT vl[3][2]; /* vertex array (last) */
86			} SPLIST; /* source polygon list */
87	greg	2.1
88	greg	2.3	extern VIEW ourview; /* our view parameters */
89			extern int hres, vres; /* our image resolution */
90			static SPLIST sphead = NULL; / our list of source polys */
91
92
93	greg	2.1	static int
94			inregion(p, cv, crit) /* check if vertex is in region */
95			FLOAT p[2];
96			double cv;
97			int crit;
98			{
99			switch (crit) {
100			case CLIP_ABOVE:
101			return(p[1] < cv);
102			case CLIP_BELOW:
103			return(p[1] >= cv);
104			case CLIP_RIGHT:
105			return(p[0] < cv);
106			case CLIP_LEFT:
107			return(p[0] >= cv);
108			}
109			return(-1);
110			}
111
112
113			static
114			clipregion(a, b, cv, crit, r) /* find intersection with boundary */
115			register FLOAT a[2], b[2];
116			double cv;
117			int crit;
118			FLOAT r[2]; /* return value */
119			{
120			switch (crit) {
121			case CLIP_ABOVE:
122			case CLIP_BELOW:
123			r[1] = cv;
124			r[0] = a[0] + (cv-a[1])/(b[1]-a[1])*(b[0]-a[0]);
125			return;
126			case CLIP_RIGHT:
127			case CLIP_LEFT:
128			r[0] = cv;
129			r[1] = a[1] + (cv-a[0])/(b[0]-a[0])*(b[1]-a[1]);
130			return;
131			}
132			}
133
134
135			static int
136			hp_clip_poly(vl, nv, cv, crit, vlo) /* clip polygon to half-plane */
137			FLOAT vl[][2];
138			int nv;
139			double cv;
140			int crit;
141			FLOAT vlo[][2]; /* return value */
142			{
143			FLOAT s, p;
144			register int j, nvo;
145
146			s = vl[nv-1];
147			nvo = 0;
148			for (j = 0; j < nv; j++) {
149			p = vl[j];
150			if (inregion(p, cv, crit)) {
151			if (!inregion(s, cv, crit))
152			clipregion(s, p, cv, crit, vlo[nvo++]);
153			vlo[nvo][0] = p[0]; vlo[nvo++][1] = p[1];
154			} else if (inregion(s, cv, crit))
155			clipregion(s, p, cv, crit, vlo[nvo++]);
156			s = p;
157			}
158			return(nvo);
159			}
160
161
162			static int
163			box_clip_poly(vl, nv, xl, xr, yb, ya, vlo) /* clip polygon to box */
164			FLOAT vl[MAXVERT][2];
165			int nv;
166			double xl, xr, yb, ya;
167			FLOAT vlo[MAXVERT][2]; /* return value */
168			{
169			FLOAT vlt[MAXVERT][2];
170			int nvt, nvo;
171
172			nvt = hp_clip_poly(vl, nv, yb, CLIP_BELOW, vlt);
173			nvo = hp_clip_poly(vlt, nvt, ya, CLIP_ABOVE, vlo);
174			nvt = hp_clip_poly(vlo, nvo, xl, CLIP_LEFT, vlt);
175			nvo = hp_clip_poly(vlt, nvt, xr, CLIP_RIGHT, vlo);
176
177			return(nvo);
178			}
179
180
181			static double
182			minw2(vl, nv, ar2) /* compute square of minimum width */
183			FLOAT vl[][2];
184			int nv;
185			double ar2;
186			{
187			double d2, w2, w2min, w2max;
188			register FLOAT p0, p1, *p2;
189			int i, j;
190			/* find minimum for all widths */
191			w2min = FHUGE;
192			p0 = vl[nv-1];
193			for (i = 0; i < nv; i++) { /* for each edge */
194			p1 = vl[i];
195			d2 = (p1[0]-p0[0])*(p1[0]-p0[0]) +
196			(p1[1]-p0[1])(p1[1]-p0[1])ar2;
197			w2max = 0.; /* find maximum for this side */
198			for (j = 1; j < nv-1; j++) {
199			p2 = vl[(i+j)%nv];
200			w2 = (p1[0]-p0[0])*(p2[1]-p0[1]) -
201			(p1[1]-p0[1])*(p2[0]-p0[0]);
202			w2 = w2w2ar2/d2; /* triangle height squared */
203			if (w2 > w2max)
204			w2max = w2;
205			}
206			if (w2max < w2min) /* global min. based on local max.'s */
207			w2min = w2max;
208			p0 = p1;
209			}
210			return(w2min);
211			}
212
213
214			static
215			convex_center(vl, nv, cv) /* compute center of convex polygon */
216			register FLOAT vl[][2];
217			int nv;
218			FLOAT cv[2]; /* return value */
219			{
220			register int i;
221			/* simple average (suboptimal) */
222			cv[0] = cv[1] = 0.;
223			for (i = 0; i < nv; i++) {
224			cv[0] += vl[i][0];
225			cv[1] += vl[i][1];
226			}
227			cv[0] /= (double)nv;
228			cv[1] /= (double)nv;
229			}
230
231
232			static double
233			poly_area(vl, nv) /* compute area of polygon */
234			register FLOAT vl[][2];
235			int nv;
236			{
237			double a;
238			FLOAT v0[2], v1[2];
239			register int i;
240
241			a = 0.;
242			v0[0] = vl[1][0] - vl[0][0];
243			v0[1] = vl[1][1] - vl[0][1];
244			for (i = 2; i < nv; i++) {
245			v1[0] = vl[i][0] - vl[0][0];
246			v1[1] = vl[i][1] - vl[0][1];
247			a += v0[0]v1[1] - v0[1]v1[0];
248			v0[0] = v1[0]; v0[1] = v1[1];
249			}
250			return(a * (a >= 0. ? .5 : -.5));
251			}
252
253
254			static int
255			convex_hull(vl, nv, vlo) /* compute polygon's convex hull */
256			FLOAT vl[][2];
257			int nv;
258			FLOAT vlo[][2]; /* return value */
259			{
260			int nvo, nvt;
261			FLOAT vlt[MAXVERT][2];
262			double voa, vta;
263			register int i, j;
264			/* start with original polygon */
265			for (i = nvo = nv; i--; ) {
266			vlo[i][0] = vl[i][0]; vlo[i][1] = vl[i][1];
267			}
268			voa = poly_area(vlo, nvo); /* compute its area */
269			for (i = 0; i < nvo; i++) { /* for each output vertex */
270			for (j = 0; j < i; j++) {
271			vlt[j][0] = vlo[j][0]; vlt[j][1] = vlo[j][1];
272			}
273			nvt = nvo - 1; /* make poly w/o vertex */
274			for (j = i; j < nvt; j++) {
275			vlt[j][0] = vlo[j+1][0]; vlt[j][1] = vlo[j+1][1];
276			}
277			vta = poly_area(vlt, nvt);
278			if (vta >= voa) { /* is simpler poly bigger? */
279			voa = vta; /* then use it */
280			for (j = nvo = nvt; j--; ) {
281			vlo[j][0] = vlt[j][0]; vlo[j][1] = vlt[j][1];
282			}
283			i--; /* next adjust */
284			}
285			}
286			return(nvo);
287			}
288
289
290	greg	2.3	static
291			spinsert(sn, vl, nv) /* insert new source polygon */
292			int sn;
293			FLOAT vl[][2];
294			int nv;
295			{
296			register SPLIST *spn;
297			register int i;
298
299			if (nv < 3)
300			return;
301			if (nv > 3)
302			spn = (SPLIST )malloc(sizeof(SPLIST)+sizeof(FLOAT)2*(nv-3));
303			else
304			spn = (SPLIST *)malloc(sizeof(SPLIST));
305			if (spn == NULL)
306			error(SYSTEM, "out of memory in spinsert");
307			spn->sn = sn;
308			for (i = spn->nv = nv; i--; ) {
309			spn->vl[i][0] = vl[i][0]; spn->vl[i][1] = vl[i][1];
310			}
311			spn->next = sphead; /* push onto global list */
312			sphead = spn;
313			}
314
315
316	greg	2.1	int
317	greg	2.3	sourcepoly(sn, sp) /* compute image polygon for source */
318	greg	2.1	int sn;
319			FLOAT sp[MAXVERT][2];
320			{
321			static char cubeord[8][6] = {{1,3,2,6,4,5},{0,4,5,7,3,2},
322			{0,1,3,7,6,4},{0,1,5,7,6,2},
323			{0,2,6,7,5,1},{0,4,6,7,3,1},
324			{0,2,3,7,5,4},{1,5,4,6,2,3}};
325			register SRCREC *s = source + sn;
326			FVECT ap, ip;
327			FLOAT pt[6][2];
328			int dir;
329			register int i, j;
330
331			if (s->sflags & (SDISTANT\|SFLAT)) {
332	greg	2.3	if (s->sflags & SDISTANT && ourview.type == VT_PAR)
333	greg	2.1	return(0); /* all or nothing case */
334			if (s->sflags & SFLAT) {
335			for (i = 0; i < 3; i++)
336	greg	2.3	ap[i] = s->sloc[i] - ourview.vp[i];
337	greg	2.1	if (DOT(ap, s->snorm) >= 0.)
338			return(0); /* source faces away */
339			}
340			for (j = 0; j < 4; j++) { /* four corners */
341			for (i = 0; i < 3; i++) {
342			ap[i] = s->sloc[i];
343			if (j==1\|j==2) ap[i] += s->ss[SU][i];
344			else ap[i] -= s->ss[SU][i];
345			if (j==2\|j==3) ap[i] += s->ss[SV][i];
346			else ap[i] -= s->ss[SV][i];
347			if (s->sflags & SDISTANT) {
348	greg	2.3	ap[i] *= 1. + ourview.vfore;
349			ap[i] += ourview.vp[i];
350	greg	2.1	}
351			}
352	greg	2.3	viewloc(ip, &ourview, ap); /* find image point */
353	greg	2.1	if (ip[2] <= 0.)
354			return(0); /* in front of view */
355			sp[j][0] = ip[0]; sp[j][1] = ip[1];
356			}
357			return(4);
358			}
359			/* identify furthest corner */
360			for (i = 0; i < 3; i++)
361	greg	2.3	ap[i] = s->sloc[i] - ourview.vp[i];
362	greg	2.1	dir = (DOT(ap,s->ss[SU])>0.) \|
363			(DOT(ap,s->ss[SV])>0.)<<1 \|
364			(DOT(ap,s->ss[SW])>0.)<<2 ;
365			/* order vertices based on this */
366			for (j = 0; j < 6; j++) {
367			for (i = 0; i < 3; i++) {
368			ap[i] = s->sloc[i];
369			if (cubeord[dir][j] & 1) ap[i] += s->ss[SU][i];
370			else ap[i] -= s->ss[SU][i];
371			if (cubeord[dir][j] & 2) ap[i] += s->ss[SV][i];
372			else ap[i] -= s->ss[SV][i];
373			if (cubeord[dir][j] & 4) ap[i] += s->ss[SW][i];
374			else ap[i] -= s->ss[SW][i];
375			}
376	greg	2.3	viewloc(ip, &ourview, ap); /* find image point */
377	greg	2.1	if (ip[2] <= 0.)
378			return(0); /* in front of view */
379			pt[j][0] = ip[0]; pt[j][1] = ip[1];
380			}
381			return(convex_hull(pt, 6, sp)); /* make sure it's convex */
382			}
383
384
385	greg	2.3	/* initialize by finding sources smaller than rad */
386			init_drawsources(rad)
387			int rad; /* source sample size */
388			{
389			FLOAT spoly[MAXVERT][2];
390			int nsv;
391			register SPLIST *sp;
392			register int i;
393			/* free old source list if one */
394			for (sp = sphead; sp != NULL; sp = sphead) {
395			sphead = sp->next;
396	greg	2.5	free((void *)sp);
397	greg	2.3	}
398			/* loop through all sources */
399			for (i = nsources; i--; ) {
400			/* compute image polygon for source */
401			if (!(nsv = sourcepoly(i, spoly)))
402			continue;
403			/* clip to image boundaries */
404			if (!(nsv = box_clip_poly(spoly, nsv, 0., 1., 0., 1., spoly)))
405			continue;
406			/* big enough for standard sampling? */
407			if (minw2(spoly, nsv, ourview.vn2/ourview.hn2) >
408			(double)rad*rad/hres/hres)
409			continue;
410			/* OK, add to our list */
411			spinsert(i, spoly, nsv);
412			}
413			}
414
415	greg	2.5	void /* add sources smaller than rad to computed subimage */
416	greg	2.3	drawsources(pic, zbf, x0, xsiz, y0, ysiz)
417	greg	2.1	COLOR pic[]; / subimage pixel value array */
418			float zbf[]; / subimage distance array (opt.) */
419			int x0, xsiz, y0, ysiz; /* origin and size of subimage */
420			{
421			FLOAT spoly[MAXVERT][2], ppoly[MAXVERT][2];
422			int nsv, npv;
423	greg	2.3	int xmin, xmax, ymin, ymax, x, y;
424	greg	2.1	FLOAT cxy[2];
425	gregl	2.4	double w;
426	greg	2.1	RAY sr;
427	greg	2.3	register SPLIST *sp;
428			register int i;
429			/* check each source in our list */
430			for (sp = sphead; sp != NULL; sp = sp->next) {
431	greg	2.1	/* clip source poly to subimage */
432	greg	2.3	nsv = box_clip_poly(sp->vl, sp->nv,
433			(double)x0/hres, (double)(x0+xsiz)/hres,
434			(double)y0/vres, (double)(y0+ysiz)/vres, spoly);
435	greg	2.1	if (!nsv)
436			continue;
437			/* find common subimage (BBox) */
438			xmin = x0 + xsiz; xmax = x0;
439			ymin = y0 + ysiz; ymax = y0;
440			for (i = 0; i < nsv; i++) {
441	greg	2.3	if ((double)xmin/hres > spoly[i][0])
442			xmin = spoly[i][0]*hres + FTINY;
443			if ((double)xmax/hres < spoly[i][0])
444			xmax = spoly[i][0]*hres - FTINY;
445			if ((double)ymin/vres > spoly[i][1])
446			ymin = spoly[i][1]*vres + FTINY;
447			if ((double)ymax/vres < spoly[i][1])
448			ymax = spoly[i][1]*vres - FTINY;
449	greg	2.1	}
450			/* evaluate each pixel in BBox */
451			for (y = ymin; y <= ymax; y++)
452			for (x = xmin; x <= xmax; x++) {
453			/* subarea for pixel */
454			npv = box_clip_poly(spoly, nsv,
455	greg	2.3	(double)x/hres, (x+1.)/hres,
456			(double)y/vres, (y+1.)/vres,
457			ppoly);
458	greg	2.1	if (!npv)
459			continue; /* no overlap */
460			convex_center(ppoly, npv, cxy);
461	greg	2.3	if ((sr.rmax = viewray(sr.rorg,sr.rdir,&ourview,
462			cxy[0],cxy[1])) < -FTINY)
463	greg	2.1	continue; /* not in view */
464	greg	2.3	if (source[sp->sn].sflags & SSPOT &&
465			spotout(&sr, source[sp->sn].sl.s))
466	greg	2.1	continue; /* outside spot */
467			rayorigin(&sr, NULL, SHADOW, 1.0);
468	greg	2.3	sr.rsrc = sp->sn;
469	greg	2.1	rayvalue(&sr); /* compute value */
470			if (bright(sr.rcol) <= FTINY)
471			continue; /* missed/blocked */
472			/* modify pixel */
473			if (zbf[y-y0] != NULL &&
474	gregl	2.4	sr.rt < 0.999*zbf[y-y0][x-x0])
475	greg	2.1	zbf[y-y0][x-x0] = sr.rt;
476	gregl	2.4	else if (!bigdiff(sr.rcol, pic[y-y0][x-x0],
477			0.001)) /* source sample */
478			setcolor(pic[y-y0][x-x0], 0., 0., 0.);
479			w = poly_area(ppoly, npv) * hres * vres;
480			scalecolor(sr.rcol, w);
481			scalecolor(pic[y-y0][x-x0], 1.-w);
482	greg	2.1	addcolor(pic[y-y0][x-x0], sr.rcol);
483			}
484			}
485			}