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 |
+ |
/* Compute average BSDF peak from current DSF's */ |
99 |
+ |
static void |
100 |
+ |
comp_bsdf_spec(void) |
101 |
+ |
{ |
102 |
+ |
const double max_hemi = 0.9; |
103 |
+ |
double peak_sum = 0; |
104 |
+ |
double rad_sum = 0; |
105 |
+ |
int n = 0; |
106 |
+ |
RBFNODE *rbf; |
107 |
+ |
FVECT sdv; |
108 |
+ |
|
109 |
+ |
if (dsf_list == NULL) { |
110 |
+ |
bsdf_spec_peak = 0; |
111 |
+ |
bsdf_spec_rad = 0; |
112 |
+ |
return; |
113 |
+ |
} |
114 |
+ |
for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) { |
115 |
+ |
sdv[0] = -rbf->invec[0]; |
116 |
+ |
sdv[1] = -rbf->invec[1]; |
117 |
+ |
sdv[2] = rbf->invec[2]*(2*(input_orient==output_orient) - 1); |
118 |
+ |
peak_sum += eval_rbfrep(rbf, sdv); |
119 |
+ |
rad_sum += est_DSFrad(rbf, sdv); |
120 |
+ |
++n; |
121 |
+ |
} |
122 |
+ |
bsdf_spec_peak = peak_sum/(double)n; |
123 |
+ |
bsdf_spec_rad = rad_sum/(double)n; |
124 |
+ |
if ((2.*M_PI)*bsdf_spec_peak*bsdf_spec_rad*bsdf_spec_rad > max_hemi) |
125 |
+ |
bsdf_spec_peak = max_hemi/((2.*M_PI)*bsdf_spec_rad*bsdf_spec_rad); |
126 |
+ |
} |
127 |
+ |
|
128 |
|
/* Create a new migration holder (sharing memory for multiprocessing) */ |
129 |
|
static MIGRATION * |
130 |
|
new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf) |
514 |
|
double best2 = M_PI*M_PI; |
515 |
|
RBFNODE *shrt_edj[2]; |
516 |
|
RBFNODE *rbf0, *rbf1; |
517 |
+ |
/* average specular peak */ |
518 |
+ |
comp_bsdf_spec(); |
519 |
|
/* add normal if needed */ |
520 |
|
check_normal_incidence(); |
521 |
|
/* check if isotropic */ |