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.31 by greg, Thu Aug 21 13:44:05 2014 UTC vs.
Revision 2.40 by greg, Tue Apr 23 14:30:36 2019 UTC

# Line 7 | Line 7 | static const char RCSid[] = "$Id$";
7   *      G. Ward
8   */
9  
10 < #ifndef _WIN32
10 > #if !defined(_WIN32) && !defined(_WIN64)
11   #include <unistd.h>
12   #include <sys/wait.h>
13   #include <sys/mman.h>
# Line 51 | Line 51 | eval_DSFsurround(const RBFNODE *rbf, const FVECT outve
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];
54 >                sum += eval_rbfrep(rbf, tvec) * COSF(tvec[2]);
55                  ++n;
56          }
57          if (n < 2)                              /* should never happen! */
# Line 63 | Line 63 | eval_DSFsurround(const RBFNODE *rbf, const FVECT outve
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];
66 >        const double    rad_epsilon = 0.01;
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 <                                                /* interpolation search */
76 <        while (outside_rad-inside_rad > rad_epsilon) {
75 >                                                /* Newton's method (sort of) */
76 >        do {
77                  double  test_rad = interp_rad;
78 <                double  DSFtest = eval_DSFsurround(rbf, outvec, test_rad);
79 <                if (DSFtarget < DSFtest) {
78 >                double  DSFtest;
79 >                if ((test_rad >= outside_rad) | (test_rad <= inside_rad))
80 >                        test_rad = .5*(inside_rad + outside_rad);
81 >                DSFtest = eval_DSFsurround(rbf, outvec, test_rad);
82 >                if (DSFtest > DSFtarget) {
83                          inside_rad = test_rad;
84                          DSFinside = DSFtest;
85                  } else {
86                          outside_rad = test_rad;
87                          DSFoutside = DSFtest;
88                  }
89 <        }
90 <        return(interp_rad);
89 >        } while (outside_rad-inside_rad > rad_epsilon);
90 >
91 >        return(.5*(inside_rad + outside_rad));
92   #undef interp_rad
93   }
94  
95 < /* Compute average BSDF peak from current DSF's */
95 > static int
96 > dbl_cmp(const void *p1, const void *p2)
97 > {
98 >        double  d1 = *(const double *)p1;
99 >        double  d2 = *(const double *)p2;
100 >
101 >        if (d1 > d2) return(1);
102 >        if (d1 < d2) return(-1);
103 >        return(0);
104 > }
105 >
106 > /* Conservative estimate of average BSDF value from current DSF's */
107   static void
108   comp_bsdf_spec(void)
109   {
110 <        double          peak_sum = 0;
110 >        double          vmod_sum = 0;
111          double          rad_sum = 0;
112          int             n = 0;
113 +        double          *cost_list = NULL;
114 +        double          max_cost = 1.;
115          RBFNODE         *rbf;
116          FVECT           sdv;
117 <
118 <        if (dsf_list == NULL) {
119 <                bsdf_spec_peak = 0;
117 >                                                /* sort by incident altitude */
118 >        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
119 >                n++;
120 >        if (n >= 10)
121 >                cost_list = (double *)malloc(sizeof(double)*n);
122 >        if (cost_list == NULL) {
123 >                bsdf_spec_val = 0;
124                  bsdf_spec_rad = 0;
125                  return;
126          }
127 +        n = 0;
128 +        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
129 +                cost_list[n++] = rbf->invec[2]*input_orient;
130 +        qsort(cost_list, n, sizeof(double), dbl_cmp);
131 +        max_cost = cost_list[(n+3)/4];          /* accept 25% nearest grazing */
132 +        free(cost_list);
133 +        n = 0;
134          for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
135 +                double  this_rad, cosfact, vest;
136 +                if (rbf->invec[2]*input_orient > max_cost)
137 +                        continue;
138                  sdv[0] = -rbf->invec[0];
139                  sdv[1] = -rbf->invec[1];
140                  sdv[2] = rbf->invec[2]*(2*(input_orient==output_orient) - 1);
141 <                peak_sum += eval_rbfrep(rbf, sdv);
142 <                rad_sum += est_DSFrad(rbf, sdv);
141 >                cosfact = COSF(sdv[2]);
142 >                this_rad = est_DSFrad(rbf, sdv);
143 >                vest = eval_rbfrep(rbf, sdv) * cosfact *
144 >                                (2.*M_PI) * this_rad*this_rad;
145 >                if (vest > rbf->vtotal)         /* don't over-estimate energy */
146 >                        vest = rbf->vtotal;
147 >                vmod_sum += vest / cosfact;     /* remove cosine factor */
148 >                rad_sum += this_rad;
149                  ++n;
150          }
114 –        bsdf_spec_peak = peak_sum/(double)n;
151          bsdf_spec_rad = rad_sum/(double)n;
152 +        bsdf_spec_val = vmod_sum/(2.*M_PI*n*bsdf_spec_rad*bsdf_spec_rad);
153   }
154  
155   /* Create a new migration holder (sharing memory for multiprocessing) */
# Line 122 | Line 159 | new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
159          size_t          memlen = sizeof(MIGRATION) +
160                                  sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
161          MIGRATION       *newmig;
162 < #ifdef _WIN32
162 > #if defined(_WIN32) || defined(_WIN64)
163          if (nprocs > 1)
164                  fprintf(stderr, "%s: warning - multiprocessing not supported\n",
165                                  progname);
# Line 153 | Line 190 | new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
190          return(mig_list = newmig);
191   }
192  
193 < #ifdef _WIN32
193 > #if defined(_WIN32) || defined(_WIN64)
194   #define await_children(n)       (void)(n)
195   #define run_subprocess()        0
196   #define end_subprocess()        (void)0
# Line 451 | Line 488 | check_normal_incidence(void)
488                  default:
489                          return;                 /* else we can interpolate */
490                  }
491 <                for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
491 >                for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
492                          const double    d = input_orient*rbf->invec[2];
493                          if (d >= 1.-2.*FTINY)
494                                  return;         /* seems we have normal */
# Line 501 | Line 538 | memerr:
538   void
539   build_mesh(void)
540   {
541 +        int             nrbfs = 0, nmigs = 0;
542          double          best2 = M_PI*M_PI;
543          RBFNODE         *shrt_edj[2];
544          RBFNODE         *rbf0, *rbf1;
545 +        const MIGRATION *ej;
546                                                  /* average specular peak */
547          comp_bsdf_spec();
548                                                  /* add normal if needed */
# Line 517 | Line 556 | build_mesh(void)
556                  return;
557          }
558          shrt_edj[0] = shrt_edj[1] = NULL;       /* start w/ shortest edge */
559 <        for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
559 >        for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next) {
560              for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
561                  double  dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
562                  if (dist2 < best2) {
# Line 525 | Line 564 | build_mesh(void)
564                          shrt_edj[1] = rbf1;
565                          best2 = dist2;
566                  }
567 +            }
568 +            ++nrbfs;
569          }
570          if (shrt_edj[0] == NULL) {
571                  fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
# Line 535 | Line 576 | build_mesh(void)
576                  mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
577          else
578                  mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
579 +                                                /* count up edges */
580 +        for (ej = mig_list; ej != NULL; ej = ej->next)
581 +                ++nmigs;
582 +        if (nmigs < nrbfs-1)                    /* did meshing fail? */
583 +                fprintf(stderr,
584 +            "%s: warning - %d incident directions but only %d interpolant(s)\n",
585 +                                progname, nrbfs, nmigs);
586                                                  /* complete migrations */
587          await_children(nchild);
588   }

Diff Legend

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