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

Comparing ray/src/cv/bsdfmesh.c (file contents):
Revision 2.27 by greg, Wed Mar 26 02:52:31 2014 UTC vs.
Revision 2.36 by greg, Sat Jan 30 17:34:00 2016 UTC

# Line 27 | Line 27 | int                    nprocs = 1;
27                                  /* number of children (-1 in child) */
28   static int              nchild = 0;
29  
30 + /* Compute average DSF value at the given radius from central vector */
31 + static double
32 + eval_DSFsurround(const RBFNODE *rbf, const FVECT outvec, const double rad)
33 + {
34 +        const int       ninc = 12;
35 +        const double    phinc = 2.*M_PI/ninc;
36 +        double          sum = 0;
37 +        int             n = 0;
38 +        FVECT           tvec;
39 +        int             i;
40 +                                                /* compute initial vector */
41 +        if (output_orient*outvec[2] >= 1.-FTINY) {
42 +                tvec[0] = tvec[2] = 0;
43 +                tvec[1] = 1;
44 +        } else {
45 +                tvec[0] = tvec[1] = 0;
46 +                tvec[2] = 1;
47 +        }
48 +        geodesic(tvec, outvec, tvec, rad, GEOD_RAD);
49 +                                                /* average surrounding DSF */
50 +        for (i = 0; i < ninc; i++) {
51 +                if (i) spinvector(tvec, tvec, outvec, phinc);
52 +                if (tvec[2] > 0 ^ output_orient > 0)
53 +                        continue;
54 +                sum += eval_rbfrep(rbf, tvec) * COSF(tvec[2]);
55 +                ++n;
56 +        }
57 +        if (n < 2)                              /* should never happen! */
58 +                return(sum);
59 +        return(sum/(double)n);
60 + }
61 +
62 + /* Estimate single-lobe radius for DSF at the given outgoing angle */
63 + static double
64 + est_DSFrad(const RBFNODE *rbf, const FVECT outvec)
65 + {
66 +        const double    rad_epsilon = 0.03;
67 +        const double    DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec) *
68 +                                                        COSF(outvec[2]);
69 +        double          inside_rad = rad_epsilon;
70 +        double          outside_rad = 0.5;
71 +        double          DSFinside = eval_DSFsurround(rbf, outvec, inside_rad);
72 +        double          DSFoutside = eval_DSFsurround(rbf, outvec, outside_rad);
73 + #define interp_rad      inside_rad + (outside_rad-inside_rad) * \
74 +                                (DSFtarget-DSFinside) / (DSFoutside-DSFinside)
75 +                                                /* Newton's method (sort of) */
76 +        do {
77 +                double  test_rad = interp_rad;
78 +                double  DSFtest;
79 +                if (test_rad >= outside_rad)
80 +                        return(test_rad);
81 +                if (test_rad <= inside_rad)
82 +                        return(test_rad*(test_rad>0));
83 +                DSFtest = eval_DSFsurround(rbf, outvec, test_rad);
84 +                if (DSFtest > DSFtarget) {
85 +                        inside_rad = test_rad;
86 +                        DSFinside = DSFtest;
87 +                } else {
88 +                        outside_rad = test_rad;
89 +                        DSFoutside = DSFtest;
90 +                }
91 +                if (DSFoutside >= DSFinside)
92 +                        return(test_rad);
93 +        } while (outside_rad-inside_rad > rad_epsilon);
94 +        return(interp_rad);
95 + #undef interp_rad
96 + }
97 +
98 + static int
99 + dbl_cmp(const void *p1, const void *p2)
100 + {
101 +        double  d1 = *(const double *)p1;
102 +        double  d2 = *(const double *)p2;
103 +
104 +        if (d1 > d2) return(1);
105 +        if (d1 < d2) return(-1);
106 +        return(0);
107 + }
108 +
109 + /* Conservative estimate of average BSDF value from current DSF's */
110 + static void
111 + comp_bsdf_spec(void)
112 + {
113 +        double          vmod_sum = 0;
114 +        double          rad_sum = 0;
115 +        int             n = 0;
116 +        double          *cost_list = NULL;
117 +        double          max_cost = 1.;
118 +        RBFNODE         *rbf;
119 +        FVECT           sdv;
120 +                                                /* sort by incident altitude */
121 +        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
122 +                n++;
123 +        if (n >= 10)
124 +                cost_list = (double *)malloc(sizeof(double)*n);
125 +        if (cost_list == NULL) {
126 +                bsdf_spec_val = 0;
127 +                bsdf_spec_rad = 0;
128 +                return;
129 +        }
130 +        n = 0;
131 +        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
132 +                cost_list[n++] = rbf->invec[2]*input_orient;
133 +        qsort(cost_list, n, sizeof(double), dbl_cmp);
134 +        max_cost = cost_list[(n+3)/4];          /* accept 25% nearest grazing */
135 +        free(cost_list);
136 +        n = 0;
137 +        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
138 +                double  this_rad, cosfact, vest;
139 +                if (rbf->invec[2]*input_orient > max_cost)
140 +                        continue;
141 +                sdv[0] = -rbf->invec[0];
142 +                sdv[1] = -rbf->invec[1];
143 +                sdv[2] = rbf->invec[2]*(2*(input_orient==output_orient) - 1);
144 +                cosfact = COSF(sdv[2]);
145 +                this_rad = est_DSFrad(rbf, sdv);
146 +                vest = eval_rbfrep(rbf, sdv) * cosfact *
147 +                                (2.*M_PI) * this_rad*this_rad;
148 +                if (vest > rbf->vtotal)         /* don't over-estimate energy */
149 +                        vest = rbf->vtotal;
150 +                vmod_sum += vest / cosfact;     /* remove cosine factor */
151 +                rad_sum += this_rad;
152 +                ++n;
153 +        }
154 +        bsdf_spec_rad = rad_sum/(double)n;
155 +        bsdf_spec_val = vmod_sum/(2.*M_PI*n*bsdf_spec_rad*bsdf_spec_rad);
156 + }
157 +
158   /* Create a new migration holder (sharing memory for multiprocessing) */
159   static MIGRATION *
160   new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
# Line 316 | Line 444 | mesh_from_edge(MIGRATION *edge)
444                                  ej1 = create_migration(tvert[0], edge->rbfv[1]);
445                          mesh_from_edge(ej0);
446                          mesh_from_edge(ej1);
447 +                        return;
448                  }
449 <        } else if (tvert[1] == NULL) {          /* grow mesh on left */
449 >        }
450 >        if (tvert[1] == NULL) {                 /* grow mesh on left */
451                  tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
452                  if (tvert[1] != NULL) {
453                          if (tvert[1]->ord > edge->rbfv[0]->ord)
# Line 394 | Line 524 | check_normal_incidence(void)
524          if (create_migration(mir_rbf, near_rbf) == NULL)
525                  exit(1);                        /* XXX should never happen! */
526          norm_vec[2] = input_orient;             /* interpolate normal dist. */
527 <        rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
527 >        rbf = e_advect_rbf(mig_list, norm_vec, 0);
528          nprocs = saved_nprocs;                  /* final clean-up */
529          free(mir_rbf);
530          free(mig_list);
# Line 414 | Line 544 | build_mesh(void)
544          double          best2 = M_PI*M_PI;
545          RBFNODE         *shrt_edj[2];
546          RBFNODE         *rbf0, *rbf1;
547 +                                                /* average specular peak */
548 +        comp_bsdf_spec();
549                                                  /* add normal if needed */
550          check_normal_incidence();
551                                                  /* check if isotropic */

Diff Legend

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