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.18 by gwlarson, Wed Nov 25 17:20:25 1998 UTC vs.
Revision 3.19 by gwlarson, Tue Dec 1 15:47:05 1998 UTC

# Line 15 | Line 15 | static char SCCSid[] = "$SunId$ SGI";
15  
16   VIEWPOINT       myeye;          /* target view position */
17  
18 + struct gclim {
19 +        HOLO    *hp;                    /* holodeck pointer */
20 +        GCOORD  gc;                     /* grid cell */
21 +        FVECT   egp;                    /* eye grid point */
22 +        double  erg[2];                 /* eye range in wall grid coords */
23 +        double  gmin[2], gmax[2];       /* grid coordinate limits */
24 + };                              /* a grid coordinate range */
25  
26 +
27 + static
28 + initeyelim(gcl, hd, gc)         /* initialize grid coordinate limits */
29 + register struct gclim   *gcl;
30 + int     hd;
31 + GCOORD  *gc;
32 + {
33 +        register FLOAT  *v;
34 +        register int    i;
35 +
36 +        gcl->hp = hdlist[hd];
37 +        copystruct(&gcl->gc, gc);
38 +        hdgrid(gcl->egp, gcl->hp, myeye.vpt);
39 +        for (i = 0; i < 2; i++) {
40 +                v = gcl->hp->wg[((gcl->gc.w>>1)+i+1)%3];
41 +                gcl->erg[i] = myeye.rng * VLEN(v);
42 +                gcl->gmin[i] = FHUGE; gcl->gmax[i] = -FHUGE;
43 +        }
44 + }
45 +
46 +
47 + static
48 + groweyelim(gcl, gp)             /* grow grid limits about eye point */
49 + register struct gclim   *gcl;
50 + FVECT   gp;
51 + {
52 +        FVECT   ab;
53 +        double  l2, d, mult, wg;
54 +        register int    i, g;
55 +
56 +        VSUB(ab, gcl->egp, gp);
57 +        l2 = DOT(ab,ab);
58 +        if (l2 <= gcl->erg[0]*gcl->erg[1]) {
59 +                gcl->gmin[0] = gcl->gmin[1] = -FHUGE;
60 +                gcl->gmax[0] = gcl->gmax[1] = FHUGE;
61 +                return;
62 +        }
63 +        mult = gp[g = gcl->gc.w>>1];
64 +        if (gcl->gc.w&1)
65 +                mult -= gcl->hp->grid[g];
66 +        if (ab[g]*ab[g] > gcl->erg[0]*gcl->erg[1])
67 +                mult /= -ab[g];
68 +        else if (fabs(ab[hdwg0[gcl->gc.w]]) > fabs(ab[hdwg1[gcl->gc.w]]))
69 +                mult = (gcl->gc.i[0] + .5 - gp[hdwg0[gcl->gc.w]]) /
70 +                                ab[hdwg0[gcl->gc.w]];
71 +        else
72 +                mult = (gcl->gc.i[1] + .5 - gp[hdwg1[gcl->gc.w]]) /
73 +                                ab[hdwg1[gcl->gc.w]];
74 +        for (i = 0; i < 2; i++) {
75 +                g = ((gcl->gc.w>>1)+i+1)%3;
76 +                wg = gp[g] + mult*ab[g];
77 +                d = mult*gcl->erg[i];
78 +                if (d < 0.) d = -d;
79 +                if (wg - d < gcl->gmin[i])
80 +                        gcl->gmin[i] = wg - d;
81 +                if (wg + d > gcl->gmax[i])
82 +                        gcl->gmax[i] = wg + d;
83 +        }
84 + }
85 +
86 +
87 + static int
88 + clipeyelim(rrng, gcl)           /* clip eye limits to grid cell */
89 + register short  rrng[2][2];
90 + register struct gclim   *gcl;
91 + {
92 +        int     incell = 1;
93 +        register int    i;
94 +
95 +        for (i = 0; i < 2; i++) {
96 +                if (gcl->gmin[i] < gcl->gc.i[i])
97 +                        gcl->gmin[i] = gcl->gc.i[i];
98 +                if (gcl->gmax[i] > gcl->gc.i[i]+1)
99 +                        gcl->gmax[i] = gcl->gc.i[i]+1;
100 +                if ((incell &= gcl->gmax[i] > gcl->gmin[i])) {
101 +                        rrng[i][0] = 256.*(gcl->gmin[i] - gcl->gc.i[i]) +
102 +                                        (1.-FTINY);
103 +                        rrng[i][1] = 256.*(gcl->gmax[i] - gcl->gc.i[i]) +
104 +                                        (1.-FTINY) - rrng[i][0];
105 +                        incell &= rrng[i][1] > 0;
106 +                }
107 +        }
108 +        return(incell);
109 + }
110 +
111 +
112   packrays(rod, p)                /* pack ray origins and directions */
113   register float  *rod;
114   register PACKET *p;
115   {
116 <        float   packdc2[RPACKSIZ];
117 <        int     iterleft = 3*p->nr + 9;
118 <        BYTE    rpos[2][2];
119 <        FVECT   ro, rd, rp1;
116 > #define gp      ro
117 > #ifdef DEBUG
118 >        double dist2sum = 0.;
119 >        FVECT   vt;
120 > #endif
121 >        int     nretries = p->nr + 2;
122 >        struct gclim    eyelim;
123 >        short   rrng0[2][2], rrng1[2][2];
124 >        int     useyelim;
125          GCOORD  gc[2];
126 <        double  d, dc2, md2, td2, dc2worst = FHUGE;
127 <        int     i;
128 <        register int    ii;
126 >        FVECT   ro, rd;
127 >        double  d;
128 >        register int    i;
129  
130          if (!hdbcoord(gc, hdlist[p->hd], p->bi))
131                  error(CONSISTENCY, "bad beam index in packrays");
132 <        td2 = myeye.rng + FTINY; td2 *= td2;
133 <        for (i = 0, md2 = 0.; i < p->nr || (md2 > td2 && iterleft--); ) {
134 <                rpos[0][0] = frandom() * 256.;
135 <                rpos[0][1] = frandom() * 256.;
136 <                rpos[1][0] = frandom() * 256.;
137 <                rpos[1][1] = frandom() * 256.;
138 <                d = hdray(ro, rd, hdlist[p->hd], gc, rpos);
139 <                if (myeye.rng > FTINY) {                /* check eyepoint */
140 <                        register int    nexti;
141 <
142 <                        VSUM(rp1, ro, rd, d);
143 <                        dc2 = dist2line(myeye.vpt, ro, rp1) / p->nr;
144 <                        if (i == p->nr) {               /* packet full */
145 <                                if (dc2 >= dc2worst)    /* quick check */
146 <                                        continue;
147 <                                nexti = 0;              /* find worst */
148 <                                for (ii = i; --ii; )
149 <                                        if (packdc2[ii] > packdc2[nexti])
150 <                                                nexti = ii;
151 <                                if (dc2 >= (dc2worst = packdc2[nexti]))
152 <                                        continue;       /* worse than worst */
153 <                                md2 -= dc2worst;
154 <                        } else
155 <                                nexti = i++;
156 <                        md2 += packdc2[nexti] = dc2;    /* new distance */
157 <                        ii = nexti;                     /* put it here */
158 <                } else
159 <                        ii = i++;
132 >        if ((useyelim = myeye.rng > FTINY)) {
133 >                initeyelim(&eyelim, p->hd, gc);
134 >                gp[gc[1].w>>1] = gc[1].w&1 ?
135 >                                hdlist[p->hd]->grid[gc[1].w>>1] : 0;
136 >                gp[hdwg0[gc[1].w]] = gc[1].i[0];
137 >                gp[hdwg1[gc[1].w]] = gc[1].i[1];
138 >                groweyelim(&eyelim, gp);
139 >                gp[hdwg0[gc[1].w]]++;
140 >                gp[hdwg1[gc[1].w]]++;
141 >                groweyelim(&eyelim, gp);
142 >                useyelim &= clipeyelim(rrng0, &eyelim);
143 >        }
144 >        for (i = 0; i < p->nr; i++) {
145 >        retry:
146 >                if (useyelim) {
147 >                        p->ra[i].r[0][0] = (int)(frandom()*rrng0[0][1])
148 >                                                + rrng0[0][0];
149 >                        p->ra[i].r[0][1] = (int)(frandom()*rrng0[1][1])
150 >                                                + rrng0[1][0];
151 >                        initeyelim(&eyelim, p->hd, gc+1);
152 >                        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);
159 >                        if (!clipeyelim(rrng1, &eyelim)) {
160 >                                useyelim &= nretries-- > 0;
161 > #ifdef DEBUG
162 >                                if (!useyelim)
163 >                                        error(WARNING, "exceeded retry limit in packrays");
164 > #endif
165 >                                goto retry;
166 >                        }
167 >                        p->ra[i].r[1][0] = (int)(frandom()*rrng1[0][1])
168 >                                                + rrng1[0][0];
169 >                        p->ra[i].r[1][1] = (int)(frandom()*rrng1[1][1])
170 >                                                + rrng1[1][0];
171 >                } else {
172 >                        p->ra[i].r[0][0] = frandom() * 256.;
173 >                        p->ra[i].r[0][1] = frandom() * 256.;
174 >                        p->ra[i].r[1][0] = frandom() * 256.;
175 >                        p->ra[i].r[1][1] = frandom() * 256.;
176 >                }
177 >                d = hdray(ro, rd, hdlist[p->hd], gc, p->ra[i].r);
178 > #ifdef DEBUG
179 >                VSUM(vt, ro, rd, d);
180 >                dist2sum += dist2line(myeye.vpt, ro, vt);
181 > #endif
182                  if (p->offset != NULL) {
183                          if (!vdef(OBSTRUCTIONS))
184                                  d *= frandom();         /* random offset */
185                          VSUM(ro, ro, rd, d);            /* advance ray */
186 <                        p->offset[ii] = d;
186 >                        p->offset[i] = d;
187                  }
188 <                p->ra[ii].r[0][0] = rpos[0][0];
189 <                p->ra[ii].r[0][1] = rpos[0][1];
190 <                p->ra[ii].r[1][0] = rpos[1][0];
191 <                p->ra[ii].r[1][1] = rpos[1][1];
72 <                VCOPY(rod+6*ii, ro);
73 <                VCOPY(rod+6*ii+3, rd);
188 >                VCOPY(rod, ro);
189 >                rod += 3;
190 >                VCOPY(rod, rd);
191 >                rod += 3;
192          }
193   #ifdef DEBUG
194 <        fprintf(stderr, "%f mean distance for target %f (%d iterations left)\n",
77 <                        sqrt(md2), myeye.rng, iterleft);
194 >        fprintf(stderr, "RMS distance = %f\n", sqrt(dist2sum/p->nr));
195   #endif
196 + #undef gp
197   }
198  
199  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines