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) * output_orient*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 |
+ |
* output_orient*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 |
+ |
/* interpolation search */ |
76 |
+ |
while (outside_rad-inside_rad > rad_epsilon) { |
77 |
+ |
double test_rad = interp_rad; |
78 |
+ |
double DSFtest = eval_DSFsurround(rbf, outvec, test_rad); |
79 |
+ |
if (DSFtarget < DSFtest) { |
80 |
+ |
inside_rad = test_rad; |
81 |
+ |
DSFinside = DSFtest; |
82 |
+ |
} else { |
83 |
+ |
outside_rad = test_rad; |
84 |
+ |
DSFoutside = DSFtest; |
85 |
+ |
} |
86 |
+ |
} |
87 |
+ |
return(interp_rad); |
88 |
+ |
#undef interp_rad |
89 |
+ |
} |
90 |
+ |
|
91 |
+ |
/* Compute average BSDF peak from current DSF's */ |
92 |
+ |
static void |
93 |
+ |
comp_bsdf_spec(void) |
94 |
+ |
{ |
95 |
+ |
double peak_sum = 0; |
96 |
+ |
double rad_sum = 0; |
97 |
+ |
int n = 0; |
98 |
+ |
RBFNODE *rbf; |
99 |
+ |
FVECT sdv; |
100 |
+ |
|
101 |
+ |
if (dsf_list == NULL) { |
102 |
+ |
bsdf_spec_peak = 0; |
103 |
+ |
bsdf_spec_crad = 0; |
104 |
+ |
return; |
105 |
+ |
} |
106 |
+ |
for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) { |
107 |
+ |
sdv[0] = -rbf->invec[0]; |
108 |
+ |
sdv[1] = -rbf->invec[1]; |
109 |
+ |
sdv[2] = rbf->invec[2]*(2*(input_orient==output_orient) - 1); |
110 |
+ |
peak_sum += eval_rbfrep(rbf, sdv); |
111 |
+ |
rad_sum += est_DSFrad(rbf, sdv); |
112 |
+ |
++n; |
113 |
+ |
} |
114 |
+ |
bsdf_spec_peak = peak_sum/(double)n; |
115 |
+ |
bsdf_spec_crad = ANG2R( rad_sum/(double)n ); |
116 |
+ |
} |
117 |
+ |
|
118 |
|
/* Create a new migration holder (sharing memory for multiprocessing) */ |
119 |
|
static MIGRATION * |
120 |
|
new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf) |
504 |
|
double best2 = M_PI*M_PI; |
505 |
|
RBFNODE *shrt_edj[2]; |
506 |
|
RBFNODE *rbf0, *rbf1; |
507 |
+ |
/* average specular peak */ |
508 |
+ |
comp_bsdf_spec(); |
509 |
|
/* add normal if needed */ |
510 |
|
check_normal_incidence(); |
511 |
|
/* check if isotropic */ |