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

Comparing ray/src/cv/bsdfinterp.c (file contents):
Revision 2.15 by greg, Wed Oct 23 03:41:39 2013 UTC vs.
Revision 2.22 by greg, Wed Sep 8 01:05:57 2021 UTC

# Line 108 | Line 108 | on_edge(const MIGRATION *ej, const FVECT ivec)
108                  return(0);
109  
110          cos_c = DOT(ej->rbfv[0]->invec, ej->rbfv[1]->invec);
111 <
112 <        return(cos_c - cos_aplusb < .001);
111 >        return(cos_c - cos_aplusb < .0002);
112   }
113  
114   /* Determine if we are inside the given triangle */
# Line 208 | Line 207 | get_interp(MIGRATION *miga[3], FVECT invec)
207          if (single_plane_incident) {            /* isotropic BSDF? */
208              RBFNODE     *rbf;                   /* find edge we're on */
209              for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
210 <                if (input_orient*rbf->invec[2] < input_orient*invec[2])
210 >                if (input_orient*rbf->invec[2] < input_orient*invec[2]-FTINY)
211                          break;
212                  if (rbf->next != NULL && input_orient*rbf->next->invec[2] <
213 <                                                        input_orient*invec[2]) {
213 >                                                input_orient*invec[2]+FTINY) {
214                      for (miga[0] = rbf->ejl; miga[0] != NULL;
215                                          miga[0] = nextedge(rbf,miga[0]))
216                          if (opp_rbf(rbf,miga[0]) == rbf->next) {
217 <                                double  nf = 1. - rbf->invec[2]*rbf->invec[2];
217 >                                double  nf = 1. -
218 >                                        rbf->next->invec[2]*rbf->next->invec[2];
219                                  if (nf > FTINY) {       /* rotate to match */
220                                          nf = sqrt((1.-invec[2]*invec[2])/nf);
221 <                                        invec[0] = nf*rbf->invec[0];
222 <                                        invec[1] = nf*rbf->invec[1];
221 >                                        invec[0] = nf*rbf->next->invec[0];
222 >                                        invec[1] = nf*rbf->next->invec[1];
223                                  }
224 <                                return(0);
224 >                                return(0);      /* rotational symmetry */
225                          }
226                      break;
227                  }
# Line 262 | Line 262 | get_interp(MIGRATION *miga[3], FVECT invec)
262          }
263   }
264  
265 < /* Advect and allocate new RBF along edge */
266 < static RBFNODE *
267 < e_advect_rbf(const MIGRATION *mig, const FVECT invec, int lobe_lim)
268 < {
269 <        double          cthresh = FTINY;
270 <        RBFNODE         *rbf;
271 <        int             n, i, j;
272 <        double          t, full_dist;
273 <                                                /* get relative position */
274 <        t = Acos(DOT(invec, mig->rbfv[0]->invec));
275 <        if (t < M_PI/grid_res) {                /* near first DSF */
276 <                n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
277 <                rbf = (RBFNODE *)malloc(n);
278 <                if (rbf == NULL)
279 <                        goto memerr;
280 <                memcpy(rbf, mig->rbfv[0], n);   /* just duplicate */
281 <                rbf->next = NULL; rbf->ejl = NULL;
282 <                return(rbf);
283 <        }
284 <        full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
285 <        if (t > full_dist-M_PI/grid_res) {      /* near second DSF */
286 <                n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
287 <                rbf = (RBFNODE *)malloc(n);
288 <                if (rbf == NULL)
289 <                        goto memerr;
290 <                memcpy(rbf, mig->rbfv[1], n);   /* just duplicate */
291 <                rbf->next = NULL; rbf->ejl = NULL;
292 <                return(rbf);
293 <        }
294 <        t /= full_dist;
295 < tryagain:
296 <        n = 0;                                  /* count migrating particles */
297 <        for (i = 0; i < mtx_nrows(mig); i++)
298 <            for (j = 0; j < mtx_ncols(mig); j++)
299 <                n += (mtx_coef(mig,i,j) > cthresh);
300 <                                                /* are we over our limit? */
301 <        if ((lobe_lim > 0) & (n > lobe_lim)) {
302 <                cthresh = cthresh*2. + 10.*FTINY;
303 <                goto tryagain;
304 <        }
305 < #ifdef DEBUG
306 <        fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n",
307 <                        mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n);
308 < #endif
309 <        rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
310 <        if (rbf == NULL)
311 <                goto memerr;
312 <        rbf->next = NULL; rbf->ejl = NULL;
313 <        VCOPY(rbf->invec, invec);
314 <        rbf->nrbf = n;
315 <        rbf->vtotal = 1.-t + t*mig->rbfv[1]->vtotal/mig->rbfv[0]->vtotal;
316 <        n = 0;                                  /* advect RBF lobes */
317 <        for (i = 0; i < mtx_nrows(mig); i++) {
318 <            const RBFVAL        *rbf0i = &mig->rbfv[0]->rbfa[i];
319 <            const float         peak0 = rbf0i->peak;
320 <            const double        rad0 = R2ANG(rbf0i->crad);
321 <            FVECT               v0;
322 <            float               mv;
323 <            ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
324 <            for (j = 0; j < mtx_ncols(mig); j++)
325 <                if ((mv = mtx_coef(mig,i,j)) > cthresh) {
326 <                        const RBFVAL    *rbf1j = &mig->rbfv[1]->rbfa[j];
327 <                        double          rad2;
328 <                        FVECT           v;
329 <                        int             pos[2];
330 <                        rad2 = R2ANG(rbf1j->crad);
331 <                        rad2 = rad0*rad0*(1.-t) + rad2*rad2*t;
332 <                        rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal *
333 <                                                rad0*rad0/rad2;
334 <                        rbf->rbfa[n].crad = ANG2R(sqrt(rad2));
335 <                        ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
336 <                        geodesic(v, v0, v, t, GEOD_REL);
337 <                        pos_from_vec(pos, v);
338 <                        rbf->rbfa[n].gx = pos[0];
339 <                        rbf->rbfa[n].gy = pos[1];
340 <                        ++n;
341 <                }
342 <        }
343 <        rbf->vtotal *= mig->rbfv[0]->vtotal;    /* turn ratio into actual */
344 <        return(rbf);
345 < memerr:
346 <        fprintf(stderr, "%s: Out of memory in e_advect_rbf()\n", progname);
347 <        exit(1);
348 <        return(NULL);   /* pro forma return */
349 < }
350 <
351 < /* Partially advect between recorded incident angles and allocate new RBF */
265 > /* Advect between recorded incident angles and allocate new RBF */
266   RBFNODE *
267   advect_rbf(const FVECT invec, int lobe_lim)
268   {
# Line 365 | Line 279 | advect_rbf(const FVECT invec, int lobe_lim)
279          VCOPY(sivec, invec);                    /* find triangle/edge */
280          sym = get_interp(miga, sivec);
281          if (sym < 0)                            /* can't interpolate? */
282 <                return(NULL);
282 >                return(def_rbf_spec(invec));
283          if (miga[1] == NULL) {                  /* advect along edge? */
284                  rbf = e_advect_rbf(miga[0], sivec, lobe_lim);
285                  if (single_plane_incident)
# Line 375 | Line 289 | advect_rbf(const FVECT invec, int lobe_lim)
289                  return(rbf);
290          }
291   #ifdef DEBUG
292 <        if (miga[0]->rbfv[0] != miga[2]->rbfv[0] |
293 <                        miga[0]->rbfv[1] != miga[1]->rbfv[0] |
294 <                        miga[1]->rbfv[1] != miga[2]->rbfv[1]) {
292 >        if ((miga[0]->rbfv[0] != miga[2]->rbfv[0]) |
293 >                        (miga[0]->rbfv[1] != miga[1]->rbfv[0]) |
294 >                        (miga[1]->rbfv[1] != miga[2]->rbfv[1])) {
295                  fprintf(stderr, "%s: Triangle vertex screw-up!\n", progname);
296                  exit(1);
297          }
# Line 428 | Line 342 | tryagain:
342              const RBFVAL        *rbf0i = &miga[0]->rbfv[0]->rbfa[i];
343              const float         w0i = rbf0i->peak;
344              const double        rad0i = R2ANG(rbf0i->crad);
345 +            C_COLOR             cc0;
346              ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
347 +            c_decodeChroma(&cc0, rbf0i->chroma);
348              for (j = 0; j < mtx_ncols(miga[0]); j++) {
349                  const float     ma = mtx_coef(miga[0],i,j);
350                  const RBFVAL    *rbf1j;
351 +                C_COLOR         ccs;
352                  double          srad2;
353                  if (ma <= cthresh)
354                          continue;
355                  rbf1j = &miga[0]->rbfv[1]->rbfa[j];
356 +                c_decodeChroma(&ccs, rbf1j->chroma);
357 +                c_cmix(&ccs, 1.-s, &cc0, s, &ccs);
358                  srad2 = R2ANG(rbf1j->crad);
359                  srad2 = (1.-s)*(1.-t)*rad0i*rad0i + s*(1.-t)*srad2*srad2;
360                  ovec_from_pos(v1, rbf1j->gx, rbf1j->gy);
# Line 453 | Line 372 | tryagain:
372                      rad2 = srad2 + t*rad2*rad2;
373                      rbf->rbfa[n].peak = w0i * ma * (mb*mbfact + mc*mcfact) *
374                                          rad0i*rad0i/rad2;
375 +                    if (rbf_colorimetry == RBCtristimulus) {
376 +                        C_COLOR cres;
377 +                        c_decodeChroma(&cres, rbf2k->chroma);
378 +                        c_cmix(&cres, 1.-t, &ccs, t, &cres);
379 +                        rbf->rbfa[n].chroma = c_encodeChroma(&cres);
380 +                    } else
381 +                        rbf->rbfa[n].chroma = c_dfchroma;
382                      rbf->rbfa[n].crad = ANG2R(sqrt(rad2));
383                      ovec_from_pos(v2, rbf2k->gx, rbf2k->gy);
384                      geodesic(v2, v1, v2, t, GEOD_REL);

Diff Legend

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