--- ray/src/cv/bsdfmesh.c 2014/08/21 19:31:23 2.32 +++ ray/src/cv/bsdfmesh.c 2017/10/06 00:23:09 2.39 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdfmesh.c,v 2.32 2014/08/21 19:31:23 greg Exp $"; +static const char RCSid[] = "$Id: bsdfmesh.c,v 2.39 2017/10/06 00:23:09 greg Exp $"; #endif /* * Create BSDF advection mesh from radial basis functions. @@ -7,7 +7,7 @@ static const char RCSid[] = "$Id: bsdfmesh.c,v 2.32 20 * G. Ward */ -#ifndef _WIN32 +#if !defined(_WIN32) && !defined(_WIN64) #include #include #include @@ -51,7 +51,7 @@ eval_DSFsurround(const RBFNODE *rbf, const FVECT outve if (i) spinvector(tvec, tvec, outvec, phinc); if (tvec[2] > 0 ^ output_orient > 0) continue; - sum += eval_rbfrep(rbf, tvec) * output_orient*tvec[2]; + sum += eval_rbfrep(rbf, tvec) * COSF(tvec[2]); ++n; } if (n < 2) /* should never happen! */ @@ -63,9 +63,9 @@ eval_DSFsurround(const RBFNODE *rbf, const FVECT outve static double est_DSFrad(const RBFNODE *rbf, const FVECT outvec) { - const double rad_epsilon = 0.03; - const double DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec) - * output_orient*outvec[2]; + const double rad_epsilon = 0.01; + const double DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec) * + COSF(outvec[2]); double inside_rad = rad_epsilon; double outside_rad = 0.5; double DSFinside = eval_DSFsurround(rbf, outvec, inside_rad); @@ -76,10 +76,8 @@ est_DSFrad(const RBFNODE *rbf, const FVECT outvec) do { double test_rad = interp_rad; double DSFtest; - if (test_rad >= outside_rad) - return(test_rad); - if (test_rad <= inside_rad) - return(test_rad*(test_rad>0)); + if ((test_rad >= outside_rad) | (test_rad <= inside_rad)) + test_rad = .5*(inside_rad + outside_rad); DSFtest = eval_DSFsurround(rbf, outvec, test_rad); if (DSFtest > DSFtarget) { inside_rad = test_rad; @@ -88,38 +86,70 @@ est_DSFrad(const RBFNODE *rbf, const FVECT outvec) outside_rad = test_rad; DSFoutside = DSFtest; } - if (DSFoutside >= DSFinside) - return(test_rad); } while (outside_rad-inside_rad > rad_epsilon); - return(interp_rad); + + return(.5*(inside_rad + outside_rad)); #undef interp_rad } -/* Compute average BSDF peak from current DSF's */ +static int +dbl_cmp(const void *p1, const void *p2) +{ + double d1 = *(const double *)p1; + double d2 = *(const double *)p2; + + if (d1 > d2) return(1); + if (d1 < d2) return(-1); + return(0); +} + +/* Conservative estimate of average BSDF value from current DSF's */ static void comp_bsdf_spec(void) { - double peak_sum = 0; + double vmod_sum = 0; double rad_sum = 0; int n = 0; + double *cost_list = NULL; + double max_cost = 1.; RBFNODE *rbf; FVECT sdv; - - if (dsf_list == NULL) { - bsdf_spec_peak = 0; + /* sort by incident altitude */ + for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) + n++; + if (n >= 10) + cost_list = (double *)malloc(sizeof(double)*n); + if (cost_list == NULL) { + bsdf_spec_val = 0; bsdf_spec_rad = 0; return; } + n = 0; + for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) + cost_list[n++] = rbf->invec[2]*input_orient; + qsort(cost_list, n, sizeof(double), dbl_cmp); + max_cost = cost_list[(n+3)/4]; /* accept 25% nearest grazing */ + free(cost_list); + n = 0; for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) { + double this_rad, cosfact, vest; + if (rbf->invec[2]*input_orient > max_cost) + continue; sdv[0] = -rbf->invec[0]; sdv[1] = -rbf->invec[1]; sdv[2] = rbf->invec[2]*(2*(input_orient==output_orient) - 1); - peak_sum += eval_rbfrep(rbf, sdv); - rad_sum += est_DSFrad(rbf, sdv); + cosfact = COSF(sdv[2]); + this_rad = est_DSFrad(rbf, sdv); + vest = eval_rbfrep(rbf, sdv) * cosfact * + (2.*M_PI) * this_rad*this_rad; + if (vest > rbf->vtotal) /* don't over-estimate energy */ + vest = rbf->vtotal; + vmod_sum += vest / cosfact; /* remove cosine factor */ + rad_sum += this_rad; ++n; } - bsdf_spec_peak = peak_sum/(double)n; bsdf_spec_rad = rad_sum/(double)n; + bsdf_spec_val = vmod_sum/(2.*M_PI*n*bsdf_spec_rad*bsdf_spec_rad); } /* Create a new migration holder (sharing memory for multiprocessing) */ @@ -129,7 +159,7 @@ new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf) size_t memlen = sizeof(MIGRATION) + sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1); MIGRATION *newmig; -#ifdef _WIN32 +#if defined(_WIN32) || defined(_WIN64) if (nprocs > 1) fprintf(stderr, "%s: warning - multiprocessing not supported\n", progname); @@ -160,7 +190,7 @@ new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf) return(mig_list = newmig); } -#ifdef _WIN32 +#if defined(_WIN32) || defined(_WIN64) #define await_children(n) (void)(n) #define run_subprocess() 0 #define end_subprocess() (void)0 @@ -458,7 +488,7 @@ check_normal_incidence(void) default: return; /* else we can interpolate */ } - for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) { + for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) { const double d = input_orient*rbf->invec[2]; if (d >= 1.-2.*FTINY) return; /* seems we have normal */