ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo2.c
(Generate patch)

Comparing ray/src/hd/rholo2.c (file contents):
Revision 3.19 by gwlarson, Tue Dec 1 15:47:05 1998 UTC vs.
Revision 3.31 by greg, Thu Apr 21 22:31:42 2022 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1998 Silicon Graphics, Inc. */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ SGI";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /*
5   * Rtrace support routines for holodeck rendering
6   */
7  
8 + #include <time.h>
9 +
10   #include "rholo.h"
11   #include "paths.h"
12   #include "random.h"
# Line 19 | Line 18 | struct gclim {
18          HOLO    *hp;                    /* holodeck pointer */
19          GCOORD  gc;                     /* grid cell */
20          FVECT   egp;                    /* eye grid point */
21 <        double  erg[2];                 /* eye range in wall grid coords */
21 >        double  erg2;                   /* mean square eye grid range */
22          double  gmin[2], gmax[2];       /* grid coordinate limits */
23   };                              /* a grid coordinate range */
24  
25 + static void initeyelim(struct gclim     *gcl, HOLO      *hp, GCOORD     *gc);
26 + static void groweyelim(struct gclim *gcl, GCOORD *gc,
27 +                double r0, double r1, int tight);
28 + static int clipeyelim(short     rrng[2][2], struct gclim        *gcl);
29  
30 < static
31 < initeyelim(gcl, hd, gc)         /* initialize grid coordinate limits */
32 < register struct gclim   *gcl;
33 < int     hd;
34 < GCOORD  *gc;
30 >
31 > static void
32 > initeyelim(             /* initialize grid coordinate limits */
33 >        struct gclim    *gcl,
34 >        HOLO    *hp,
35 >        GCOORD  *gc
36 > )
37   {
38 <        register FLOAT  *v;
39 <        register int    i;
38 >        RREAL   *v;
39 >        int     i;
40  
41 <        gcl->hp = hdlist[hd];
42 <        copystruct(&gcl->gc, gc);
43 <        hdgrid(gcl->egp, gcl->hp, myeye.vpt);
44 <        for (i = 0; i < 2; i++) {
45 <                v = gcl->hp->wg[((gcl->gc.w>>1)+i+1)%3];
46 <                gcl->erg[i] = myeye.rng * VLEN(v);
42 <                gcl->gmin[i] = FHUGE; gcl->gmax[i] = -FHUGE;
41 >        if (hp != NULL) {
42 >                hdgrid(gcl->egp, gcl->hp = hp, myeye.vpt);
43 >                gcl->erg2 = 0;
44 >                for (i = 0, v = hp->wg[0]; i < 3; i++, v += 3)
45 >                        gcl->erg2 += DOT(v,v);
46 >                gcl->erg2 *= (1./3.) * myeye.rng*myeye.rng;
47          }
48 +        if (gc != NULL)
49 +                gcl->gc = *gc;
50 +        gcl->gmin[0] = gcl->gmin[1] = FHUGE;
51 +        gcl->gmax[0] = gcl->gmax[1] = -FHUGE;
52   }
53  
54  
55 < static
56 < groweyelim(gcl, gp)             /* grow grid limits about eye point */
57 < register struct gclim   *gcl;
58 < FVECT   gp;
55 > static void
56 > groweyelim(     /* grow grid limits about eye point */
57 >        struct gclim    *gcl,
58 >        GCOORD  *gc,
59 >        double  r0,
60 >        double  r1,
61 >        int     tight
62 > )
63   {
64 <        FVECT   ab;
65 <        double  l2, d, mult, wg;
66 <        register int    i, g;
67 <
64 >        FVECT   gp, ab;
65 >        double  ab2, od, cfact;
66 >        double  sqcoef[3], ctcoef[3], licoef[3], cnst;
67 >        int     gw, gi[2];
68 >        double  wallpos, a, b, c, d, e, f;
69 >        double  root[2], yex;
70 >        int     n, i, j, nex;
71 >                                                /* point/view cone */
72 >        i = gc->w>>1;
73 >        gp[i] = gc->w&1 ? gcl->hp->grid[i] : 0;
74 >        gp[hdwg0[gc->w]] = gc->i[0] + r0;
75 >        gp[hdwg1[gc->w]] = gc->i[1] + r1;
76          VSUB(ab, gcl->egp, gp);
77 <        l2 = DOT(ab,ab);
78 <        if (l2 <= gcl->erg[0]*gcl->erg[1]) {
77 >        ab2 = DOT(ab, ab);
78 >        gw = gcl->gc.w>>1;
79 >        if ((i==gw ? ab[gw]*ab[gw] : ab2)  <= gcl->erg2 + FTINY) {
80                  gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
81                  gcl->gmax[0] = gcl->gmax[1] = FHUGE;
82 <                return;
82 >                return;                 /* too close (to wall) */
83          }
84 <        mult = gp[g = gcl->gc.w>>1];
85 <        if (gcl->gc.w&1)
86 <                mult -= gcl->hp->grid[g];
87 <        if (ab[g]*ab[g] > gcl->erg[0]*gcl->erg[1])
88 <                mult /= -ab[g];
89 <        else if (fabs(ab[hdwg0[gcl->gc.w]]) > fabs(ab[hdwg1[gcl->gc.w]]))
90 <                mult = (gcl->gc.i[0] + .5 - gp[hdwg0[gcl->gc.w]]) /
91 <                                ab[hdwg0[gcl->gc.w]];
92 <        else
93 <                mult = (gcl->gc.i[1] + .5 - gp[hdwg1[gcl->gc.w]]) /
94 <                                ab[hdwg1[gcl->gc.w]];
84 >        ab2 = 1./ab2;                           /* 1/norm2(ab) */
85 >        od = DOT(gp, ab);                       /* origin dot direction */
86 >        cfact = 1./(1. - ab2*gcl->erg2);        /* tan^2 + 1 of cone angle */
87 >        for (i = 0; i < 3; i++) {               /* compute cone equation */
88 >                sqcoef[i] = ab[i]*ab[i]*cfact*ab2 - 1.;
89 >                ctcoef[i] = 2.*ab[i]*ab[(i+1)%3]*cfact*ab2;
90 >                licoef[i] = 2.*(gp[i] - ab[i]*cfact*od*ab2);
91 >        }
92 >        cnst = cfact*od*od*ab2 - DOT(gp,gp);
93 >        /*
94 >         * CONE:        sqcoef[0]*x*x + sqcoef[1]*y*y + sqcoef[2]*z*z
95 >         *              + ctcoef[0]*x*y + ctcoef[1]*y*z + ctcoef[2]*z*x
96 >         *              + licoef[0]*x + licoef[1]*y + licoef[2]*z + cnst == 0
97 >         */
98 >                                /* equation for conic section in plane */
99 >        gi[0] = hdwg0[gcl->gc.w];
100 >        gi[1] = hdwg1[gcl->gc.w];
101 >        wallpos = gcl->gc.w&1 ? gcl->hp->grid[gw] : 0;
102 >        a = sqcoef[gi[0]];                                      /* x2 */
103 >        b = ctcoef[gi[0]];                                      /* xy */
104 >        c = sqcoef[gi[1]];                                      /* y2 */
105 >        d = ctcoef[gw]*wallpos + licoef[gi[0]];                 /* x */
106 >        e = ctcoef[gi[1]]*wallpos + licoef[gi[1]];              /* y */
107 >        f = wallpos*(wallpos*sqcoef[gw] + licoef[gw]) + cnst;
108          for (i = 0; i < 2; i++) {
109 <                g = ((gcl->gc.w>>1)+i+1)%3;
110 <                wg = gp[g] + mult*ab[g];
111 <                d = mult*gcl->erg[i];
112 <                if (d < 0.) d = -d;
113 <                if (wg - d < gcl->gmin[i])
114 <                        gcl->gmin[i] = wg - d;
115 <                if (wg + d > gcl->gmax[i])
116 <                        gcl->gmax[i] = wg + d;
109 >                if (i) {                /* swap x and y coefficients */
110 >                        double  t;
111 >                        t = a; a = c; c = t;
112 >                        t = d; d = e; e = t;
113 >                }
114 >                nex = 0;                /* check global extrema */
115 >                n = quadratic(root, a*(4.*a*c-b*b), 2.*a*(2.*c*d-b*e),
116 >                                d*(c*d-b*e) + f*b*b);
117 >                while (n-- > 0) {
118 >                        if (gc->w>>1 == gi[i] &&
119 >                                        (gc->w&1) ^ (root[n] < gp[gc->w>>1])) {
120 >                                if (gc->w&1)
121 >                                        gcl->gmin[i] = -FHUGE;
122 >                                else
123 >                                        gcl->gmax[i] = FHUGE;
124 >                                nex++;
125 >                                continue;               /* hyperbolic */
126 >                        }
127 >                        if (tight) {
128 >                                yex = (-2.*a*root[n] - d)/b;
129 >                                if (yex < gcl->gc.i[1-i] ||
130 >                                                yex > gcl->gc.i[1-i]+1)
131 >                                        continue;       /* outside cell */
132 >                        }
133 >                        if (root[n] < gcl->gmin[i])
134 >                                gcl->gmin[i] = root[n];
135 >                        if (root[n] > gcl->gmax[i])
136 >                                gcl->gmax[i] = root[n];
137 >                        nex++;
138 >                }
139 >                                        /* check local extrema */
140 >                for (j = nex < 2 ? 2 : 0; j--; ) {
141 >                        yex = gcl->gc.i[1-i] + j;
142 >                        n = quadratic(root, a, b*yex+d, yex*(yex*c+e)+f);
143 >                        while (n-- > 0) {
144 >                                if (gc->w>>1 == gi[i] &&
145 >                                        (gc->w&1) ^ (root[n] < gp[gc->w>>1]))
146 >                                        continue;
147 >                                if (root[n] < gcl->gmin[i])
148 >                                        gcl->gmin[i] = root[n];
149 >                                if (root[n] > gcl->gmax[i])
150 >                                        gcl->gmax[i] = root[n];
151 >                        }
152 >                }
153          }
154   }
155  
156  
157   static int
158 < clipeyelim(rrng, gcl)           /* clip eye limits to grid cell */
159 < register short  rrng[2][2];
160 < register struct gclim   *gcl;
158 > clipeyelim(             /* clip eye limits to grid cell */
159 >        short   rrng[2][2],
160 >        struct gclim    *gcl
161 > )
162   {
163          int     incell = 1;
164 <        register int    i;
164 >        int     i;
165  
166          for (i = 0; i < 2; i++) {
167                  if (gcl->gmin[i] < gcl->gc.i[i])
168                          gcl->gmin[i] = gcl->gc.i[i];
169                  if (gcl->gmax[i] > gcl->gc.i[i]+1)
170                          gcl->gmax[i] = gcl->gc.i[i]+1;
171 <                if ((incell &= gcl->gmax[i] > gcl->gmin[i])) {
171 >                if (gcl->gmax[i] > gcl->gmin[i]) {
172                          rrng[i][0] = 256.*(gcl->gmin[i] - gcl->gc.i[i]) +
173                                          (1.-FTINY);
174                          rrng[i][1] = 256.*(gcl->gmax[i] - gcl->gc.i[i]) +
175                                          (1.-FTINY) - rrng[i][0];
176 <                        incell &= rrng[i][1] > 0;
177 <                }
176 >                } else
177 >                        rrng[i][0] = rrng[i][1] = 0;
178 >                incell &= rrng[i][1] > 0;
179          }
180          return(incell);
181   }
182  
183  
184 < packrays(rod, p)                /* pack ray origins and directions */
185 < register float  *rod;
186 < register PACKET *p;
184 > void
185 > packrays(               /* pack ray origins and directions */
186 >        float   *rod,
187 >        PACKET  *p
188 > )
189   {
190 < #define gp      ro
191 < #ifdef DEBUG
118 <        double dist2sum = 0.;
190 > #if 0
191 >        double  dist2sum = 0.;
192          FVECT   vt;
193   #endif
194          int     nretries = p->nr + 2;
# Line 125 | Line 198 | register PACKET        *p;
198          GCOORD  gc[2];
199          FVECT   ro, rd;
200          double  d;
201 <        register int    i;
201 >        int     i;
202  
203          if (!hdbcoord(gc, hdlist[p->hd], p->bi))
204                  error(CONSISTENCY, "bad beam index in packrays");
205          if ((useyelim = myeye.rng > FTINY)) {
206 <                initeyelim(&eyelim, p->hd, gc);
207 <                gp[gc[1].w>>1] = gc[1].w&1 ?
208 <                                hdlist[p->hd]->grid[gc[1].w>>1] : 0;
209 <                gp[hdwg0[gc[1].w]] = gc[1].i[0];
210 <                gp[hdwg1[gc[1].w]] = gc[1].i[1];
211 <                groweyelim(&eyelim, gp);
212 <                gp[hdwg0[gc[1].w]]++;
213 <                gp[hdwg1[gc[1].w]]++;
141 <                groweyelim(&eyelim, gp);
142 <                useyelim &= clipeyelim(rrng0, &eyelim);
206 >                initeyelim(&eyelim, hdlist[p->hd], gc);
207 >                groweyelim(&eyelim, gc+1, 0., 0., 0);
208 >                groweyelim(&eyelim, gc+1, 1., 1., 0);
209 >                useyelim = clipeyelim(rrng0, &eyelim);
210 > #ifdef DEBUG
211 >                if (!useyelim)
212 >                        error(WARNING, "no eye overlap in packrays");
213 > #endif
214          }
215          for (i = 0; i < p->nr; i++) {
216          retry:
217                  if (useyelim) {
218 <                        p->ra[i].r[0][0] = (int)(frandom()*rrng0[0][1])
219 <                                                + rrng0[0][0];
220 <                        p->ra[i].r[0][1] = (int)(frandom()*rrng0[1][1])
221 <                                                + rrng0[1][0];
222 <                        initeyelim(&eyelim, p->hd, gc+1);
223 <                        gp[gc[0].w>>1] = gc[0].w&1 ?
153 <                                        hdlist[p->hd]->grid[gc[0].w>>1] : 0;
154 <                        gp[hdwg0[gc[0].w]] = gc[0].i[0] +
155 <                                        (1./256.)*(p->ra[i].r[0][0]+.5);
156 <                        gp[hdwg1[gc[0].w]] = gc[0].i[1] +
157 <                                        (1./256.)*(p->ra[i].r[0][1]+.5);
158 <                        groweyelim(&eyelim, gp);
218 >                        initeyelim(&eyelim, NULL, gc+1);
219 >                        p->ra[i].r[0][0] = irandom(rrng0[0][1]) + rrng0[0][0];
220 >                        p->ra[i].r[0][1] = irandom(rrng0[1][1]) + rrng0[1][0];
221 >                        groweyelim(&eyelim, gc,
222 >                                        (1./256.)*(p->ra[i].r[0][0]+.5),
223 >                                        (1./256.)*(p->ra[i].r[0][1]+.5), 1);
224                          if (!clipeyelim(rrng1, &eyelim)) {
225 <                                useyelim &= nretries-- > 0;
225 >                                useyelim = nretries-- > 0;
226   #ifdef DEBUG
227                                  if (!useyelim)
228 <                                        error(WARNING, "exceeded retry limit in packrays");
228 >                                        error(WARNING,
229 >                                        "exceeded retry limit in packrays");
230   #endif
231                                  goto retry;
232                          }
233 <                        p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1])
234 <                                                + rrng1[0][0];
169 <                        p->ra[i].r[1][1] = (int)(frandom()*rrng1[1][1])
170 <                                                + rrng1[1][0];
233 >                        p->ra[i].r[1][0] = irandom(rrng1[0][1]) + rrng1[0][0];
234 >                        p->ra[i].r[1][1] = irandom(rrng1[1][1]) + rrng1[1][0];
235                  } else {
236 <                        p->ra[i].r[0][0] = frandom() * 256.;
237 <                        p->ra[i].r[0][1] = frandom() * 256.;
238 <                        p->ra[i].r[1][0] = frandom() * 256.;
239 <                        p->ra[i].r[1][1] = frandom() * 256.;
236 >                        p->ra[i].r[0][0] = random() & 0xff;
237 >                        p->ra[i].r[0][1] = random() & 0xff;
238 >                        p->ra[i].r[1][0] = random() & 0xff;
239 >                        p->ra[i].r[1][1] = random() & 0xff;
240                  }
241                  d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
242 < #ifdef DEBUG
242 > #if 0
243                  VSUM(vt, ro, rd, d);
244                  dist2sum += dist2line(myeye.vpt, ro, vt);
245   #endif
# Line 190 | Line 254 | register PACKET        *p;
254                  VCOPY(rod, rd);
255                  rod += 3;
256          }
257 < #ifdef DEBUG
258 <        fprintf(stderr, "RMS distance = %f\n", sqrt(dist2sum/p->nr));
257 > #if 0
258 >        fprintf(stderr, "%f RMS (%d retries)\t", sqrt(dist2sum/p->nr),
259 >                        p->nr + 2 - nretries);
260   #endif
196 – #undef gp
261   }
262  
263  
264 < donerays(p, rvl)                /* encode finished ray computations */
265 < register PACKET *p;
266 < register float  *rvl;
264 > void
265 > donerays(               /* encode finished ray computations */
266 >        PACKET  *p,
267 >        float   *rvl
268 > )
269   {
270          double  d;
271 <        register int    i;
271 >        int     i;
272  
273          for (i = 0; i < p->nr; i++) {
274                  setcolr(p->ra[i].v, rvl[0], rvl[1], rvl[2]);
# Line 217 | Line 283 | register float *rvl;
283  
284  
285   int
286 < done_rtrace()                   /* clean up and close rtrace calculation */
286 > done_rtrace(void)                       /* clean up and close rtrace calculation */
287   {
288          int     status;
289                                          /* already closed? */
290          if (!nprocs)
291 <                return;
291 >                return(0);
292                                          /* flush beam queue */
293          done_packets(flush_queue());
294                                          /* sync holodeck */
# Line 238 | Line 304 | done_rtrace()                  /* clean up and close rtrace calculati
304   }
305  
306  
307 < new_rtrace()                    /* restart rtrace calculation */
307 > void
308 > new_rtrace(void)                        /* restart rtrace calculation */
309   {
310          char    combuf[128];
311  
# Line 262 | Line 329 | new_rtrace()                   /* restart rtrace calculation */
329   }
330  
331  
332 < getradfile()                    /* run rad and get needed variables */
332 > int
333 > getradfile(void)                        /* run rad and get needed variables */
334   {
335          static short    mvar[] = {OCTREE,EYESEP,-1};
336          static char     tf1[] = TEMPLATE;
337          char    tf2[64];
338          char    combuf[256];
339 <        char    *pippt;
340 <        register int    i;
341 <        register char   *cp;
339 >        char    *pippt = NULL;
340 >        int     i;
341 >        char    *cp;
342                                          /* check if rad file specified */
343          if (!vdef(RIF))
344                  return(0);
# Line 293 | Line 361 | getradfile()                   /* run rad and get needed variables */
361                          pippt = NULL;
362                  }
363          if (pippt != NULL)
364 <                strcpy(pippt, "> /dev/null");   /* nothing to match */
364 >                strcpy(pippt, "> " NULL_DEVICE);        /* nothing to match */
365          else
366                  sprintf(cp, ")[ \t]*=' > %s", tf2);
367   #ifdef DEBUG
# Line 304 | Line 372 | getradfile()                   /* run rad and get needed variables */
372                  loadvars(tf2);                  /* load variables */
373                  unlink(tf2);
374          }
375 <        rtargc += wordfile(rtargv+rtargc, tf1); /* get rtrace options */
375 >                                                /* get rtrace options */
376 >        rtargc += wordfile(rtargv+rtargc, MAXRTARGC-rtargc, tf1);
377          unlink(tf1);                    /* clean up */
378          return(1);
379   }
380  
381  
382 < report(t)                       /* report progress so far */
383 < time_t  t;
382 > void
383 > report(                 /* report progress so far */
384 >        time_t  t
385 > )
386   {
387          static time_t   seconds2go = 1000000;
388  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines