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.1 by greg, Fri Oct 19 04:14:29 2012 UTC vs.
Revision 2.10 by greg, Thu Sep 26 14:57:18 2013 UTC

# Line 23 | Line 23 | int                    nprocs = 1;
23                                  /* number of children (-1 in child) */
24   static int              nchild = 0;
25  
26 < /* Compute (and allocate) migration price matrix for optimization */
27 < static float *
28 < price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
29 < {
30 <        float   *pmtx = (float *)malloc(sizeof(float) *
31 <                                        from_rbf->nrbf * to_rbf->nrbf);
32 <        FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
33 <        int     i, j;
26 > typedef struct {
27 >        int             nrows, ncols;   /* array size (matches migration) */
28 >        float           *price;         /* migration prices */
29 >        short           *sord;          /* sort for each row, low to high */
30 >        float           *prow;          /* current price row */
31 > } PRICEMAT;                     /* sorted pricing matrix */
32  
33 <        if ((pmtx == NULL) | (vto == NULL)) {
34 <                fprintf(stderr, "%s: Out of memory in migration_costs()\n",
37 <                                progname);
38 <                exit(1);
39 <        }
40 <        for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
41 <                ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
33 > #define pricerow(p,i)   ((p)->price + (i)*(p)->ncols)
34 > #define psortrow(p,i)   ((p)->sord + (i)*(p)->ncols)
35  
43 –        for (i = from_rbf->nrbf; i--; ) {
44 –            const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
45 –            FVECT               vfrom;
46 –            ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
47 –            for (j = to_rbf->nrbf; j--; )
48 –                pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
49 –                                fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
50 –        }
51 –        free(vto);
52 –        return(pmtx);
53 – }
54 –
55 – /* Comparison routine needed for sorting price row */
56 – static const float      *price_arr;
57 – static int
58 – msrt_cmp(const void *p1, const void *p2)
59 – {
60 –        float   c1 = price_arr[*(const int *)p1];
61 –        float   c2 = price_arr[*(const int *)p2];
62 –
63 –        if (c1 > c2) return(1);
64 –        if (c1 < c2) return(-1);
65 –        return(0);
66 – }
67 –
68 – /* Compute minimum (optimistic) cost for moving the given source material */
69 – static double
70 – min_cost(double amt2move, const double *avail, const float *price, int n)
71 – {
72 –        static int      *price_sort = NULL;
73 –        static int      n_alloc = 0;
74 –        double          total_cost = 0;
75 –        int             i;
76 –
77 –        if (amt2move <= FTINY)                  /* pre-emptive check */
78 –                return(0.);
79 –        if (n > n_alloc) {                      /* (re)allocate sort array */
80 –                if (n_alloc) free(price_sort);
81 –                price_sort = (int *)malloc(sizeof(int)*n);
82 –                if (price_sort == NULL) {
83 –                        fprintf(stderr, "%s: Out of memory in min_cost()\n",
84 –                                        progname);
85 –                        exit(1);
86 –                }
87 –                n_alloc = n;
88 –        }
89 –        for (i = n; i--; )
90 –                price_sort[i] = i;
91 –        price_arr = price;
92 –        qsort(price_sort, n, sizeof(int), &msrt_cmp);
93 –                                                /* move cheapest first */
94 –        for (i = 0; i < n && amt2move > FTINY; i++) {
95 –                int     d = price_sort[i];
96 –                double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
97 –
98 –                total_cost += amt * price[d];
99 –                amt2move -= amt;
100 –        }
101 –        return(total_cost);
102 – }
103 –
104 – /* Take a step in migration by choosing optimal bucket to transfer */
105 – static double
106 – migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
107 – {
108 –        const double    maxamt = .1;
109 –        const double    minamt = maxamt*.0001;
110 –        static double   *src_cost = NULL;
111 –        static int      n_alloc = 0;
112 –        struct {
113 –                int     s, d;   /* source and destination */
114 –                double  price;  /* price estimate per amount moved */
115 –                double  amt;    /* amount we can move */
116 –        } cur, best;
117 –        int             i;
118 –
119 –        if (mtx_nrows(mig) > n_alloc) {         /* allocate cost array */
120 –                if (n_alloc)
121 –                        free(src_cost);
122 –                src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
123 –                if (src_cost == NULL) {
124 –                        fprintf(stderr, "%s: Out of memory in migration_step()\n",
125 –                                        progname);
126 –                        exit(1);
127 –                }
128 –                n_alloc = mtx_nrows(mig);
129 –        }
130 –        for (i = mtx_nrows(mig); i--; )         /* starting costs for diff. */
131 –                src_cost[i] = min_cost(src_rem[i], dst_rem,
132 –                                        pmtx+i*mtx_ncols(mig), mtx_ncols(mig));
133 –
134 –                                                /* find best source & dest. */
135 –        best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
136 –        for (cur.s = mtx_nrows(mig); cur.s--; ) {
137 –            const float *price = pmtx + cur.s*mtx_ncols(mig);
138 –            double      cost_others = 0;
139 –            if (src_rem[cur.s] < minamt)
140 –                    continue;
141 –            cur.d = -1;                         /* examine cheapest dest. */
142 –            for (i = mtx_ncols(mig); i--; )
143 –                if (dst_rem[i] > minamt &&
144 –                                (cur.d < 0 || price[i] < price[cur.d]))
145 –                        cur.d = i;
146 –            if (cur.d < 0)
147 –                    return(.0);
148 –            if ((cur.price = price[cur.d]) >= best.price)
149 –                    continue;                   /* no point checking further */
150 –            cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
151 –                                src_rem[cur.s] : dst_rem[cur.d];
152 –            if (cur.amt > maxamt) cur.amt = maxamt;
153 –            dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
154 –            for (i = mtx_nrows(mig); i--; )
155 –                if (i != cur.s)
156 –                        cost_others += min_cost(src_rem[i], dst_rem,
157 –                                                price, mtx_ncols(mig))
158 –                                        - src_cost[i];
159 –            dst_rem[cur.d] += cur.amt;          /* undo trial move */
160 –            cur.price += cost_others/cur.amt;   /* adjust effective price */
161 –            if (cur.price < best.price)         /* are we better than best? */
162 –                    best = cur;
163 –        }
164 –        if ((best.s < 0) | (best.d < 0))
165 –                return(.0);
166 –                                                /* make the actual move */
167 –        mig->mtx[mtx_ndx(mig,best.s,best.d)] += best.amt;
168 –        src_rem[best.s] -= best.amt;
169 –        dst_rem[best.d] -= best.amt;
170 –        return(best.amt);
171 – }
172 –
173 – #ifdef DEBUG
174 – static char *
175 – thetaphi(const FVECT v)
176 – {
177 –        static char     buf[128];
178 –        double          theta, phi;
179 –
180 –        theta = 180./M_PI*acos(v[2]);
181 –        phi = 180./M_PI*atan2(v[1],v[0]);
182 –        sprintf(buf, "(%.0f,%.0f)", theta, phi);
183 –
184 –        return(buf);
185 – }
186 – #endif
187 –
36   /* Create a new migration holder (sharing memory for multiprocessing) */
37   static MIGRATION *
38   new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
# Line 272 | Line 120 | run_subprocess(void)
120                  if (pid < 0) {
121                          fprintf(stderr, "%s: cannot fork subprocess\n",
122                                          progname);
123 +                        await_children(nchild);
124                          exit(1);
125                  }
126                  ++nchild;                       /* subprocess started */
# Line 286 | Line 135 | run_subprocess(void)
135  
136   #endif  /* ! _WIN32 */
137  
138 + /* Comparison routine needed for sorting price row */
139 + static int
140 + msrt_cmp(void *b, const void *p1, const void *p2)
141 + {
142 +        PRICEMAT        *pm = (PRICEMAT *)b;
143 +        float           c1 = pm->prow[*(const short *)p1];
144 +        float           c2 = pm->prow[*(const short *)p2];
145 +
146 +        if (c1 > c2) return(1);
147 +        if (c1 < c2) return(-1);
148 +        return(0);
149 + }
150 +
151 + /* Compute (and allocate) migration price matrix for optimization */
152 + static void
153 + price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
154 + {
155 +        FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
156 +        int     i, j;
157 +
158 +        pm->nrows = from_rbf->nrbf;
159 +        pm->ncols = to_rbf->nrbf;
160 +        pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
161 +        pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
162 +        
163 +        if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
164 +                fprintf(stderr, "%s: Out of memory in migration_costs()\n",
165 +                                progname);
166 +                exit(1);
167 +        }
168 +        for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
169 +                ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
170 +
171 +        for (i = from_rbf->nrbf; i--; ) {
172 +            const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
173 +            FVECT               vfrom;
174 +            short               *srow;
175 +            ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
176 +            pm->prow = pricerow(pm,i);
177 +            srow = psortrow(pm,i);
178 +            for (j = to_rbf->nrbf; j--; ) {
179 +                double          dprod = DOT(vfrom, vto[j]);
180 +                pm->prow[j] = ((dprod >= 1.) ? .0 : acos(dprod)) +
181 +                                fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
182 +                srow[j] = j;
183 +            }
184 +            qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
185 +        }
186 +        free(vto);
187 + }
188 +
189 + /* Free price matrix */
190 + static void
191 + free_routes(PRICEMAT *pm)
192 + {
193 +        free(pm->price); pm->price = NULL;
194 +        free(pm->sord); pm->sord = NULL;
195 + }
196 +
197 + /* Compute minimum (optimistic) cost for moving the given source material */
198 + static double
199 + min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
200 + {
201 +        double          total_cost = 0;
202 +        int             j;
203 +
204 +        if (amt2move <= FTINY)                  /* pre-emptive check */
205 +                return(.0);
206 +                                                /* move cheapest first */
207 +        for (j = 0; j < pm->ncols && amt2move > FTINY; j++) {
208 +                int     d = psortrow(pm,s)[j];
209 +                double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
210 +
211 +                total_cost += amt * pricerow(pm,s)[d];
212 +                amt2move -= amt;
213 +        }
214 +        return(total_cost);
215 + }
216 +
217 + /* Take a step in migration by choosing optimal bucket to transfer */
218 + static double
219 + migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const PRICEMAT *pm)
220 + {
221 +        const double    maxamt = 1./(double)pm->ncols;
222 +        const double    minamt = maxamt*5e-6;
223 +        double          *src_cost;
224 +        struct {
225 +                int     s, d;   /* source and destination */
226 +                double  price;  /* price estimate per amount moved */
227 +                double  amt;    /* amount we can move */
228 +        } cur, best;
229 +        int             i;
230 +                                                /* allocate cost array */
231 +        src_cost = (double *)malloc(sizeof(double)*pm->nrows);
232 +        if (src_cost == NULL) {
233 +                fprintf(stderr, "%s: Out of memory in migration_step()\n",
234 +                                progname);
235 +                exit(1);
236 +        }
237 +        for (i = pm->nrows; i--; )              /* starting costs for diff. */
238 +                src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
239 +
240 +                                                /* find best source & dest. */
241 +        best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
242 +        for (cur.s = pm->nrows; cur.s--; ) {
243 +            double      cost_others = 0;
244 +
245 +            if (src_rem[cur.s] <= minamt)
246 +                    continue;
247 +                                                /* examine cheapest dest. */
248 +            for (i = 0; i < pm->ncols; i++)
249 +                if (dst_rem[ cur.d = psortrow(pm,cur.s)[i] ] > minamt)
250 +                        break;
251 +            if (i >= pm->ncols)
252 +                break;
253 +            if ((cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price)
254 +                continue;                       /* no point checking further */
255 +            cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
256 +                                src_rem[cur.s] : dst_rem[cur.d];
257 +            if (cur.amt > maxamt) cur.amt = maxamt;
258 +            dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
259 +            for (i = pm->nrows; i--; )
260 +                if (i != cur.s)
261 +                        cost_others += min_cost(src_rem[i], dst_rem, pm, i)
262 +                                        - src_cost[i];
263 +            dst_rem[cur.d] += cur.amt;          /* undo trial move */
264 +            cur.price += cost_others/cur.amt;   /* adjust effective price */
265 +            if (cur.price < best.price)         /* are we better than best? */
266 +                    best = cur;
267 +        }
268 +        free(src_cost);                         /* finish up */
269 +
270 +        if ((best.s < 0) | (best.d < 0))        /* nothing left to move? */
271 +                return(.0);
272 +                                                /* else make the actual move */
273 +        mtx_coef(mig,best.s,best.d) += best.amt;
274 +        src_rem[best.s] -= best.amt;
275 +        dst_rem[best.d] -= best.amt;
276 +        return(best.amt);
277 + }
278 +
279   /* Compute and insert migration along directed edge (may fork child) */
280   static MIGRATION *
281   create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
282   {
283 <        const double    end_thresh = 0.1/(from_rbf->nrbf*to_rbf->nrbf);
284 <        const double    check_thresh = 0.01;
295 <        const double    rel_thresh = 5e-6;
296 <        float           *pmtx;
283 >        const double    end_thresh = 5e-6;
284 >        PRICEMAT        pmtx;
285          MIGRATION       *newmig;
286          double          *src_rem, *dst_rem;
287          double          total_rem = 1., move_amt;
288 <        int             i;
288 >        int             i, j;
289                                                  /* check if exists already */
290          for (newmig = from_rbf->ejl; newmig != NULL;
291                          newmig = nextedge(from_rbf,newmig))
292                  if (newmig->rbfv[1] == to_rbf)
293                          return(NULL);
294                                                  /* else allocate */
295 + #ifdef DEBUG
296 +        fprintf(stderr, "Building path from (theta,phi) (%.0f,%.0f) ",
297 +                        get_theta180(from_rbf->invec),
298 +                        get_phi360(from_rbf->invec));
299 +        fprintf(stderr, "to (%.0f,%.0f) with %d x %d matrix\n",
300 +                        get_theta180(to_rbf->invec),
301 +                        get_phi360(to_rbf->invec),
302 +                        from_rbf->nrbf, to_rbf->nrbf);
303 + #endif
304          newmig = new_migration(from_rbf, to_rbf);
305          if (run_subprocess())
306                  return(newmig);                 /* child continues */
307 <        pmtx = price_routes(from_rbf, to_rbf);
307 >        price_routes(&pmtx, from_rbf, to_rbf);
308          src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
309          dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
310          if ((src_rem == NULL) | (dst_rem == NULL)) {
# Line 315 | Line 312 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
312                                  progname);
313                  exit(1);
314          }
318 – #ifdef DEBUG
319 –        fprintf(stderr, "Building path from (theta,phi) %s ",
320 –                        thetaphi(from_rbf->invec));
321 –        fprintf(stderr, "to %s", thetaphi(to_rbf->invec));
322 –        /* if (nchild) */ fputc('\n', stderr);
323 – #endif
315                                                  /* starting quantities */
316          memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
317          for (i = from_rbf->nrbf; i--; )
318                  src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
319 <        for (i = to_rbf->nrbf; i--; )
320 <                dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
319 >        for (j = to_rbf->nrbf; j--; )
320 >                dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
321 >
322          do {                                    /* move a bit at a time */
323 <                move_amt = migration_step(newmig, src_rem, dst_rem, pmtx);
323 >                move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
324                  total_rem -= move_amt;
325 < #ifdef DEBUG
326 <                if (!nchild)
335 <                        /* fputc('.', stderr); */
336 <                        fprintf(stderr, "%.9f remaining...\r", total_rem);
337 < #endif
338 <        } while (total_rem > end_thresh && (total_rem > check_thresh) |
339 <                                        (move_amt > rel_thresh*total_rem));
340 < #ifdef DEBUG
341 <        if (!nchild) fputs("\ndone.\n", stderr);
342 <        else fprintf(stderr, "finished with %.9f remaining\n", total_rem);
343 < #endif
325 >        } while ((total_rem > end_thresh) & (move_amt > 0));
326 >
327          for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
328 <            float       nf = rbf_volume(&from_rbf->rbfa[i]);
346 <            int         j;
328 >            double      nf = rbf_volume(&from_rbf->rbfa[i]);
329              if (nf <= FTINY) continue;
330              nf = from_rbf->vtotal / nf;
331              for (j = to_rbf->nrbf; j--; )
332 <                newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
332 >                mtx_coef(newmig,i,j) *= nf;     /* row now sums to 1.0 */
333          }
334          end_subprocess();                       /* exit here if subprocess */
335 <        free(pmtx);                             /* free working arrays */
335 >        free_routes(&pmtx);                     /* free working arrays */
336          free(src_rem);
337          free(dst_rem);
338          return(newmig);

Diff Legend

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