ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/rt/ambcomp.c
(Generate patch)

Comparing ray/src/rt/ambcomp.c (file contents):
Revision 2.93 by greg, Tue Apr 16 23:32:20 2024 UTC vs.
Revision 2.100 by greg, Mon Apr 28 19:30:01 2025 UTC

# Line 24 | Line 24 | static const char      RCSid[] = "$Id$";
24   #ifndef MINADIV
25   #define MINADIV         7       /* minimum # divisions in each dimension */
26   #endif
27 + #ifndef MINSDIST
28 + #define MINSDIST        0.25    /* def. min. spacing = 1/4th division */
29 + #endif
30  
31   typedef struct {
32          FVECT   p;              /* intersection point */
# Line 62 | Line 65 | ambcollision(                          /* proposed direciton collides? */
65   {
66          double  cos_thresh;
67          int     ii, jj;
68 <                                        /* min. spacing = 1/4th division */
69 <        cos_thresh = (PI/4.)/(double)hp->ns;
68 >
69 >        cos_thresh = (PI*MINSDIST)/(double)hp->ns;
70          cos_thresh = 1. - .5*cos_thresh*cos_thresh;
71                                          /* check existing neighbors */
72          for (ii = i-1; ii <= i+1; ii++) {
# Line 89 | Line 92 | ambcollision(                          /* proposed direciton collides? */
92   }
93  
94  
95 + #define XLOTSIZ         251             /* size of used car lot */
96 + #define CFIRST          0               /* first corner */
97 + #define COTHER          (CFIRST+4)      /* non-corner sample */
98 + #define CMAXTARGET      (int)(XLOTSIZ*MINSDIST/(1-MINSDIST))
99 + #define CXCOPY(d,s)     (excharr[d][0]=excharr[s][0], excharr[d][1]=excharr[s][1])
100 +
101   static int
102 + psample_class(double ss[2])             /* classify patch sample */
103 + {
104 +        if (ss[0] < MINSDIST) {
105 +                if (ss[1] < MINSDIST)
106 +                        return(CFIRST);
107 +                if (ss[1] > 1.-MINSDIST)
108 +                        return(CFIRST+2);
109 +        } else if (ss[0] > 1.-MINSDIST) {
110 +                if (ss[1] < MINSDIST)
111 +                        return(CFIRST+1);
112 +                if (ss[1] > 1.-MINSDIST)
113 +                        return(CFIRST+3);
114 +        }
115 +        return(COTHER);                 /* not in a corner */
116 + }
117 +
118 + static void
119 + trade_patchsamp(double ss[2])           /* trade in problem patch position */
120 + {
121 +        static float    excharr[XLOTSIZ][2];
122 +        static short    gterm[COTHER+1];
123 +        double          srep[2];
124 +        int             sclass, rclass;
125 +        int             x;
126 +                                        /* reset on corner overload */
127 +        if (gterm[COTHER-1] >= (CMAXTARGET+XLOTSIZ)/2)
128 +                memset(gterm, 0, sizeof(gterm));
129 +                                        /* (re-)initialize? */
130 +        while (gterm[COTHER] < XLOTSIZ) {
131 +                excharr[gterm[COTHER]][0] = frandom();
132 +                excharr[gterm[COTHER]][1] = frandom();
133 +                ++gterm[COTHER];
134 +        }                               /* get trade-in candidate... */
135 +        sclass = psample_class(ss);     /* submitted corner or not? */
136 +        switch (sclass) {
137 +        case COTHER:                    /* trade mid-edge with corner/any */
138 +                x = irandom( gterm[COTHER-1] > CMAXTARGET
139 +                                ? gterm[COTHER-1] : XLOTSIZ );
140 +                break;
141 +        case CFIRST:                    /* kick out of first corner */
142 +                x = gterm[CFIRST] + irandom(XLOTSIZ - gterm[CFIRST]);
143 +                break;
144 +        default:                        /* kick out of 2nd-4th corner */
145 +                x = irandom(XLOTSIZ - (gterm[sclass] - gterm[sclass-1]));
146 +                x += (x >= gterm[sclass-1])*(gterm[sclass] - gterm[sclass-1]);
147 +                break;
148 +        }
149 +        srep[0] = excharr[x][0];        /* save selected replacement (result) */
150 +        srep[1] = excharr[x][1];
151 +                                        /* identify replacement class */
152 +        for (rclass = CFIRST; rclass < COTHER; rclass++)
153 +                if (x < gterm[rclass])
154 +                        break;          /* repark to keep classes grouped */
155 +        while (rclass > sclass) {       /* replacement group after submitted? */
156 +                CXCOPY(x, gterm[rclass-1]);
157 +                x = gterm[--rclass]++;
158 +        }
159 +        while (rclass < sclass) {       /* replacement group before submitted? */
160 +                --gterm[rclass];
161 +                CXCOPY(x, gterm[rclass]);
162 +                x = gterm[rclass++];
163 +        }
164 +        excharr[x][0] = ss[0];          /* complete the trade-in */
165 +        excharr[x][1] = ss[1];
166 +        ss[0] = srep[0];
167 +        ss[1] = srep[1];
168 + }
169 +
170 + #undef CXCOPY
171 + #undef XLOTSIZ
172 + #undef COTHER
173 + #undef CFIRST
174 +
175 +
176 + static int
177   ambsample(                              /* initial ambient division sample */
178          AMBHEMI *hp,
179          int     i,
# Line 100 | Line 184 | ambsample(                             /* initial ambient division sample */
184          AMBSAMP *ap = &ambsam(hp,i,j);
185          RAY     ar;
186          int     hlist[3], ii;
187 +        double  ss[2];
188          RREAL   spt[2];
189          double  zd;
190                                          /* generate hemispherical sample */
# Line 115 | Line 200 | ambsample(                             /* initial ambient division sample */
200                  scalescolor(ar.rcoef, 1./AVGREFL);
201          }
202          hlist[0] = hp->rp->rno;
203 <        hlist[1] = j;
204 <        hlist[2] = i;
205 <        multisamp(spt, 2, urand(ilhash(hlist,3)+n));
206 < resample:
207 <        square2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns);
203 >        hlist[1] = AI(hp,i,j);
204 >        hlist[2] = samplendx;
205 >        multisamp(ss, 2, urand(ilhash(hlist,3)+n));
206 > patch_redo:
207 >        square2disk(spt, (j+ss[1])/hp->ns, (i+ss[0])/hp->ns);
208          zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]);
209          for (ii = 3; ii--; )
210                  ar.rdir[ii] =   spt[0]*hp->ux[ii] +
# Line 127 | Line 212 | resample:
212                                  zd*hp->onrm[ii];
213          checknorm(ar.rdir);
214                                          /* avoid coincident samples */
215 <        if (!n && ambcollision(hp, i, j, ar.rdir)) {
216 <                spt[0] = frandom(); spt[1] = frandom();
217 <                goto resample;          /* reject this sample */
215 >        if (!n & (hp->ns >= 4) && ambcollision(hp, i, j, ar.rdir)) {
216 >                trade_patchsamp(ss);
217 >                goto patch_redo;
218          }
219          dimlist[ndims++] = AI(hp,i,j) + 90171;
220          rayvalue(&ar);                  /* evaluate ray */
# Line 163 | Line 248 | static float *
248   getambdiffs(AMBHEMI *hp)
249   {
250          const double    normf = 1./(pbright(hp->acoef) + FTINY);
251 <        float   *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float));
252 <        float   *ep, *earr2;
251 >        float   *earr = (float *)calloc(2*hp->ns*hp->ns, sizeof(float));
252 >        float   *ep;
253          AMBSAMP *ap;
254          double  b, b1, d2;
255          int     i, j;
# Line 172 | Line 257 | getambdiffs(AMBHEMI *hp)
257          if (earr == NULL)               /* out of memory? */
258                  return(NULL);
259                                          /* sum squared neighbor diffs */
260 <        for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++)
260 >        ap = hp->sa;
261 >        ep = earr + hp->ns*hp->ns;      /* original estimates to scratch */
262 >        for (i = 0; i < hp->ns; i++)
263              for (j = 0; j < hp->ns; j++, ap++, ep++) {
264                  b = pbright(ap[0].v);
265                  if (i) {                /* from above */
# Line 198 | Line 285 | getambdiffs(AMBHEMI *hp)
285                  ep[-hp->ns-1] += d2;
286              }
287                                          /* correct for number of neighbors */
288 <        earr[0] *= 6./3.;
289 <        earr[hp->ns-1] *= 6./3.;
290 <        earr[(hp->ns-1)*hp->ns] *= 6./3.;
291 <        earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 6./3.;
288 >        ep = earr + hp->ns*hp->ns;
289 >        ep[0] *= 6./3.;
290 >        ep[hp->ns-1] *= 6./3.;
291 >        ep[(hp->ns-1)*hp->ns] *= 6./3.;
292 >        ep[(hp->ns-1)*hp->ns + hp->ns-1] *= 6./3.;
293          for (i = 1; i < hp->ns-1; i++) {
294 <                earr[i*hp->ns] *= 6./5.;
295 <                earr[i*hp->ns + hp->ns-1] *= 6./5.;
294 >                ep[i*hp->ns] *= 6./5.;
295 >                ep[i*hp->ns + hp->ns-1] *= 6./5.;
296          }
297          for (j = 1; j < hp->ns-1; j++) {
298 <                earr[j] *= 6./5.;
299 <                earr[(hp->ns-1)*hp->ns + j] *= 6./5.;
298 >                ep[j] *= 6./5.;
299 >                ep[(hp->ns-1)*hp->ns + j] *= 6./5.;
300          }
301 <                                        /* preen map to avoid cliffs */
214 <        earr2 = (float *)malloc(hp->ns*hp->ns*sizeof(float));
215 <        if (earr2 == NULL)
216 <                return(earr);
217 <        memcpy(earr2, earr, hp->ns*hp->ns*sizeof(float));
301 >                                        /* blur final map to reduce bias */
302          for (i = 0; i < hp->ns-1; i++) {
303 <            float  *ep2 = earr2 + i*hp->ns;
303 >            float  *ep2;
304              ep = earr + i*hp->ns;
305 <            for (j = 0; j < hp->ns-1; j++, ep2++, ep++) {
306 <                if (ep2[1] < .5*ep2[0]) {
307 <                        ep[0] -= .125*ep2[0];
308 <                        ep[1] += .125*ep2[0];
309 <                } else if (ep2[1] > 2.*ep2[0]) {
226 <                        ep[1] -= .125*ep2[1];
227 <                        ep[0] += .125*ep2[1];
228 <                }
229 <                if (ep2[hp->ns] < .5*ep2[0]) {
230 <                        ep[0] -= .125*ep2[0];
231 <                        ep[hp->ns] += .125*ep2[0];
232 <                } else if (ep2[hp->ns] > 2.*ep2[0]) {
233 <                        ep[hp->ns] -= .125*ep2[hp->ns];
234 <                        ep[0] += .125*ep2[hp->ns];
235 <                }
305 >            ep2 = ep + hp->ns*hp->ns;
306 >            for (j = 0; j < hp->ns-1; j++, ep++, ep2++) {
307 >                ep[0] += .5*ep2[0] + .125*(ep2[1] + ep2[hp->ns]);
308 >                ep[1] += .125*ep2[0];
309 >                ep[hp->ns] += .125*ep2[0];
310              }
311          }
238 –        free(earr2);
312          return(earr);
313   }
314  
# Line 287 | Line 360 | samp_hemi(                             /* sample indirect hemisphere */
360                                          /* set number of divisions */
361          if (backside) wt = -wt;
362          if (ambacc <= FTINY &&
363 <                        wt > (d *= 0.8*r->rweight/(ambdiv*minweight)))
363 >                        wt > (d *= 0.8*r->rweight/(ambdiv*minweight + 1e-20)))
364                  wt = d;                 /* avoid ray termination */
365          n = sqrt(ambdiv * wt) + 0.5;
366          i = 1 + (MINADIV-1)*(ambacc > FTINY);
# Line 336 | Line 409 | samp_hemi(                             /* sample indirect hemisphere */
409          if (hp->sampOK <= MINADIV*MINADIV)
410                  return(hp);             /* don't bother super-sampling */
411          n = ambssamp*wt + 0.5;
412 <        if (n > 8) {                    /* perform super-sampling? */
412 >        if (n >= 4*hp->ns) {            /* perform super-sampling? */
413                  ambsupersamp(hp, n);
414                  copyscolor(rcol, hp->acol);
415          }

Diff Legend

– Removed lines
+ Added lines
< Changed lines (old)
> Changed lines (new)