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.14 by gwlarson, Mon Nov 23 17:50:26 1998 UTC vs.
Revision 3.18 by gwlarson, Wed Nov 25 17:20:25 1998 UTC

# Line 13 | Line 13 | static char SCCSid[] = "$SunId$ SGI";
13   #include "random.h"
14  
15  
16 + VIEWPOINT       myeye;          /* target view position */
17 +
18 +
19   packrays(rod, p)                /* pack ray origins and directions */
20 < float   *rod;
20 > register float  *rod;
21   register PACKET *p;
22   {
23 <        static int      nmh = 0;
24 <        static int      *mhtab;
25 <        FVECT   ro, rd;
26 <        register BEAM   *bp;
23 >        float   packdc2[RPACKSIZ];
24 >        int     iterleft = 3*p->nr + 9;
25 >        BYTE    rpos[2][2];
26 >        FVECT   ro, rd, rp1;
27          GCOORD  gc[2];
28 <        int     ila[4], offset;
29 <        double  d, sl[4];
30 <        register int    i, j, k;
28 >        double  d, dc2, md2, td2, dc2worst = FHUGE;
29 >        int     i;
30 >        register int    ii;
31  
32          if (!hdbcoord(gc, hdlist[p->hd], p->bi))
33                  error(CONSISTENCY, "bad beam index in packrays");
34 <                                                        /* uniqueness hash */
35 <        if ((bp = hdgetbeam(hdlist[p->hd], p->bi)) != NULL) {
36 <                if (2*bp->nrm > nmh) {
37 <                        if (nmh) free((char *)mhtab);
38 <                        nmh = 2*bp->nrm + 1;
39 <                        mhtab = (int *)malloc(nmh*sizeof(int));
40 <                        if (mhtab == NULL)
41 <                                error(SYSTEM, "out of memory in packrays");
42 <                }
43 <                for (k = nmh; k--; )
44 <                        mhtab[k] = -1;
45 <                for (i = bp->nrm; i--; ) {
46 <                        ila[0] = hdbray(bp)[i].r[0][0];
47 <                        ila[1] = hdbray(bp)[i].r[0][1];
48 <                        ila[2] = hdbray(bp)[i].r[1][0];
49 <                        ila[3] = hdbray(bp)[i].r[1][1];
50 <                        for (k = ilhash(ila,4); mhtab[k%nmh] >= 0; k++)
51 <                                ;
52 <                        mhtab[k%nmh] = i;
53 <                }
54 <        }
55 <                                                        /* init each ray */
56 <        ila[0] = p->hd; ila[1] = p->bi;
57 <        offset = ilhash(ila,2) + p->nc;
58 <        for (i = 0; i < p->nr; i++) {
59 <                do {                                    /* next unique ray */
60 <                        multisamp(sl, 4, urand(offset+i));
61 <                        p->ra[i].r[0][0] = ila[0] = sl[0] * 256.;
59 <                        p->ra[i].r[0][1] = ila[1] = sl[1] * 256.;
60 <                        p->ra[i].r[1][0] = ila[2] = sl[2] * 256.;
61 <                        p->ra[i].r[1][1] = ila[3] = sl[3] * 256.;
62 <                        if (bp == NULL)
63 <                                break;
64 <                        for (k = ilhash(ila,4); (j = mhtab[k%nmh]) >= 0; k++)
65 <                                if (hdbray(bp)[j].r[0][0] ==
66 <                                                        p->ra[i].r[0][0] &&
67 <                                                hdbray(bp)[j].r[0][1] ==
68 <                                                        p->ra[i].r[0][1] &&
69 <                                                hdbray(bp)[j].r[1][0] ==
70 <                                                        p->ra[i].r[1][0] &&
71 <                                                hdbray(bp)[j].r[1][1] ==
72 <                                                        p->ra[i].r[1][1]) {
73 <                                        offset += bp->nrm - j;
74 <                                        break;
75 <                                }
76 <                } while (j >= 0);
77 <                d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
78 <                if (!vdef(OBSTRUCTIONS))
79 <                        d *= frandom();                 /* random offset */
34 >        td2 = myeye.rng + FTINY; td2 *= td2;
35 >        for (i = 0, md2 = 0.; i < p->nr || (md2 > td2 && iterleft--); ) {
36 >                rpos[0][0] = frandom() * 256.;
37 >                rpos[0][1] = frandom() * 256.;
38 >                rpos[1][0] = frandom() * 256.;
39 >                rpos[1][1] = frandom() * 256.;
40 >                d = hdray(ro, rd, hdlist[p->hd], gc, rpos);
41 >                if (myeye.rng > FTINY) {                /* check eyepoint */
42 >                        register int    nexti;
43 >
44 >                        VSUM(rp1, ro, rd, d);
45 >                        dc2 = dist2line(myeye.vpt, ro, rp1) / p->nr;
46 >                        if (i == p->nr) {               /* packet full */
47 >                                if (dc2 >= dc2worst)    /* quick check */
48 >                                        continue;
49 >                                nexti = 0;              /* find worst */
50 >                                for (ii = i; --ii; )
51 >                                        if (packdc2[ii] > packdc2[nexti])
52 >                                                nexti = ii;
53 >                                if (dc2 >= (dc2worst = packdc2[nexti]))
54 >                                        continue;       /* worse than worst */
55 >                                md2 -= dc2worst;
56 >                        } else
57 >                                nexti = i++;
58 >                        md2 += packdc2[nexti] = dc2;    /* new distance */
59 >                        ii = nexti;                     /* put it here */
60 >                } else
61 >                        ii = i++;
62                  if (p->offset != NULL) {
63 +                        if (!vdef(OBSTRUCTIONS))
64 +                                d *= frandom();         /* random offset */
65                          VSUM(ro, ro, rd, d);            /* advance ray */
66 <                        p->offset[i] = d;
66 >                        p->offset[ii] = d;
67                  }
68 <                VCOPY(rod, ro);
69 <                rod += 3;
70 <                VCOPY(rod, rd);
71 <                rod += 3;
68 >                p->ra[ii].r[0][0] = rpos[0][0];
69 >                p->ra[ii].r[0][1] = rpos[0][1];
70 >                p->ra[ii].r[1][0] = rpos[1][0];
71 >                p->ra[ii].r[1][1] = rpos[1][1];
72 >                VCOPY(rod+6*ii, ro);
73 >                VCOPY(rod+6*ii+3, rd);
74          }
75 + #ifdef DEBUG
76 +        fprintf(stderr, "%f mean distance for target %f (%d iterations left)\n",
77 +                        sqrt(md2), myeye.rng, iterleft);
78 + #endif
79   }
80  
81  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines