ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/pabopto2xml.c
Revision: 2.10
Committed: Tue Sep 18 02:50:13 2012 UTC (11 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.9: +455 -53 lines
Log Message: 
Further additions in preparation for actual interpolation

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.10 static const char RCSid[] = "$Id: pabopto2xml.c,v 2.9 2012/09/06 00:07:43 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * Convert PAB-Opto measurements to XML format using tensor tree representation
6     * Employs Bonneel et al. Earth Mover's Distance interpolant.
7     *
8     * G.Ward
9     */
10    
11     #define _USE_MATH_DEFINES
12     #include <stdio.h>
13     #include <stdlib.h>
14     #include <string.h>
15     #include <ctype.h>
16     #include <math.h>
17     #include "bsdf.h"
18    
19 greg 2.10 #define DEBUG 1
20    
21 greg 2.1 #ifndef GRIDRES
22 greg 2.10 #define GRIDRES 200 /* grid resolution per side */
23 greg 2.1 #endif
24    
25 greg 2.3 #define RSCA 2.7 /* radius scaling factor (empirical) */
26 greg 2.2
27 greg 2.6 /* convert to/from coded radians */
28     #define ANG2R(r) (int)((r)*((1<<16)/M_PI))
29 greg 2.2 #define R2ANG(c) (((c)+.5)*(M_PI/(1<<16)))
30 greg 2.1
31     typedef struct {
32 greg 2.5 float vsum; /* DSF sum */
33 greg 2.1 unsigned short nval; /* number of values in sum */
34 greg 2.2 unsigned short crad; /* radius (coded angle) */
35 greg 2.1 } GRIDVAL; /* grid value */
36    
37     typedef struct {
38 greg 2.5 float peak; /* lobe value at peak */
39 greg 2.2 unsigned short crad; /* radius (coded angle) */
40 greg 2.1 unsigned char gx, gy; /* grid position */
41     } RBFVAL; /* radial basis function value */
42    
43 greg 2.7 struct s_rbfnode; /* forward declaration of RBF struct */
44    
45     typedef struct s_migration {
46     struct s_migration *next; /* next in global edge list */
47     struct s_rbfnode *rbfv[2]; /* from,to vertex */
48     struct s_migration *enxt[2]; /* next from,to sibling */
49     float mtx[1]; /* matrix (extends struct) */
50     } MIGRATION; /* migration link (winged edge structure) */
51    
52     typedef struct s_rbfnode {
53     struct s_rbfnode *next; /* next in global RBF list */
54     MIGRATION *ejl; /* edge list for this vertex */
55 greg 2.1 FVECT invec; /* incident vector direction */
56 greg 2.8 double vtotal; /* volume for normalization */
57 greg 2.1 int nrbf; /* number of RBFs */
58     RBFVAL rbfa[1]; /* RBF array (extends struct) */
59 greg 2.10 } RBFNODE; /* RBF representation of DSF @ 1 incidence */
60 greg 2.1
61     /* our loaded grid for this incident angle */
62 greg 2.10 static double theta_in_deg, phi_in_deg;
63     static GRIDVAL dsf_grid[GRIDRES][GRIDRES];
64    
65     /* all incident angles in-plane so far? */
66     static int single_plane_incident = -1;
67    
68     /* input/output orientations */
69     static int input_orient = 0;
70     static int output_orient = 0;
71 greg 2.1
72 greg 2.5 /* processed incident DSF measurements */
73 greg 2.10 static RBFNODE *dsf_list = NULL;
74 greg 2.7
75 greg 2.8 /* RBF-linking matrices (edges) */
76 greg 2.7 static MIGRATION *mig_list = NULL;
77    
78 greg 2.10 /* migration edges drawn in raster fashion */
79     static MIGRATION *mig_grid[GRIDRES][GRIDRES];
80    
81 greg 2.8 #define mtx_nrows(m) ((m)->rbfv[0]->nrbf)
82     #define mtx_ncols(m) ((m)->rbfv[1]->nrbf)
83     #define mtx_ndx(m,i,j) ((i)*mtx_ncols(m) + (j))
84     #define is_src(rbf,m) ((rbf) == (m)->rbfv[0])
85     #define is_dest(rbf,m) ((rbf) == (m)->rbfv[1])
86     #define nextedge(rbf,m) (m)->enxt[is_dest(rbf,m)]
87 greg 2.10 #define opp_rbf(rbf,m) (m)->rbfv[is_src(rbf,m)]
88    
89     #define round(v) (int)((v) + .5 - ((v) < -.5))
90 greg 2.8
91     /* Compute volume associated with Gaussian lobe */
92     static double
93     rbf_volume(const RBFVAL *rbfp)
94     {
95     double rad = R2ANG(rbfp->crad);
96    
97     return((2.*M_PI) * rbfp->peak * rad*rad);
98     }
99 greg 2.1
100 greg 2.3 /* Compute outgoing vector from grid position */
101     static void
102 greg 2.10 ovec_from_pos(FVECT vec, int xpos, int ypos)
103 greg 2.1 {
104 greg 2.3 double uv[2];
105     double r2;
106    
107     SDsquare2disk(uv, (1./GRIDRES)*(xpos+.5), (1./GRIDRES)*(ypos+.5));
108     /* uniform hemispherical projection */
109     r2 = uv[0]*uv[0] + uv[1]*uv[1];
110     vec[0] = vec[1] = sqrt(2. - r2);
111     vec[0] *= uv[0];
112     vec[1] *= uv[1];
113 greg 2.10 vec[2] = output_orient*(1. - r2);
114 greg 2.1 }
115    
116 greg 2.10 /* Compute grid position from normalized input/output vector */
117 greg 2.1 static void
118     pos_from_vec(int pos[2], const FVECT vec)
119     {
120     double sq[2]; /* uniform hemispherical projection */
121 greg 2.10 double norm = 1./sqrt(1. + fabs(vec[2]));
122 greg 2.1
123     SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
124    
125     pos[0] = (int)(sq[0]*GRIDRES);
126     pos[1] = (int)(sq[1]*GRIDRES);
127     }
128    
129 greg 2.5 /* Evaluate RBF for DSF at the given normalized outgoing direction */
130 greg 2.1 static double
131 greg 2.10 eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
132 greg 2.1 {
133     double res = .0;
134     const RBFVAL *rbfp;
135     FVECT odir;
136     double sig2;
137     int n;
138    
139     rbfp = rp->rbfa;
140     for (n = rp->nrbf; n--; rbfp++) {
141 greg 2.10 ovec_from_pos(odir, rbfp->gx, rbfp->gy);
142 greg 2.2 sig2 = R2ANG(rbfp->crad);
143     sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
144 greg 2.1 if (sig2 > -19.)
145 greg 2.5 res += rbfp->peak * exp(sig2);
146 greg 2.1 }
147     return(res);
148     }
149    
150 greg 2.10 /* Insert a new directional scattering function in our global list */
151     static void
152     insert_dsf(RBFNODE *newrbf)
153     {
154     RBFNODE *rbf, *rbf_last;
155    
156     /* keep in ascending theta order */
157     for (rbf_last = NULL, rbf = dsf_list;
158     single_plane_incident & (rbf != NULL);
159     rbf_last = rbf, rbf = rbf->next)
160     if (input_orient*rbf->invec[2] < input_orient*newrbf->invec[2])
161     break;
162     if (rbf_last == NULL) {
163     newrbf->next = dsf_list;
164     dsf_list = newrbf;
165     return;
166     }
167     newrbf->next = rbf;
168     rbf_last->next = newrbf;
169     }
170    
171 greg 2.3 /* Count up filled nodes and build RBF representation from current grid */
172 greg 2.10 static RBFNODE *
173 greg 2.3 make_rbfrep(void)
174     {
175 greg 2.6 int niter = 16;
176     double lastVar, thisVar = 100.;
177 greg 2.3 int nn;
178 greg 2.10 RBFNODE *newnode;
179 greg 2.3 int i, j;
180    
181     nn = 0; /* count selected bins */
182     for (i = 0; i < GRIDRES; i++)
183     for (j = 0; j < GRIDRES; j++)
184 greg 2.6 nn += dsf_grid[i][j].nval;
185 greg 2.3 /* allocate RBF array */
186 greg 2.10 newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
187 greg 2.3 if (newnode == NULL) {
188 greg 2.8 fputs("Out of memory in make_rbfrep()\n", stderr);
189 greg 2.3 exit(1);
190     }
191     newnode->next = NULL;
192 greg 2.7 newnode->ejl = NULL;
193 greg 2.3 newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
194     newnode->invec[0] = cos(M_PI/180.*phi_in_deg)*newnode->invec[2];
195     newnode->invec[1] = sin(M_PI/180.*phi_in_deg)*newnode->invec[2];
196 greg 2.10 newnode->invec[2] = input_orient*sqrt(1. - newnode->invec[2]*newnode->invec[2]);
197 greg 2.8 newnode->vtotal = 0;
198 greg 2.3 newnode->nrbf = nn;
199     nn = 0; /* fill RBF array */
200     for (i = 0; i < GRIDRES; i++)
201     for (j = 0; j < GRIDRES; j++)
202 greg 2.5 if (dsf_grid[i][j].nval) {
203 greg 2.6 newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
204 greg 2.5 newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
205 greg 2.3 newnode->rbfa[nn].gx = i;
206     newnode->rbfa[nn].gy = j;
207     ++nn;
208     }
209 greg 2.6 /* iterate to improve interpolation accuracy */
210     do {
211 greg 2.4 double dsum = .0, dsum2 = .0;
212 greg 2.3 nn = 0;
213     for (i = 0; i < GRIDRES; i++)
214     for (j = 0; j < GRIDRES; j++)
215 greg 2.5 if (dsf_grid[i][j].nval) {
216 greg 2.3 FVECT odir;
217 greg 2.6 double corr;
218 greg 2.10 ovec_from_pos(odir, i, j);
219 greg 2.6 newnode->rbfa[nn++].peak *= corr =
220 greg 2.5 dsf_grid[i][j].vsum /
221 greg 2.3 eval_rbfrep(newnode, odir);
222 greg 2.4 dsum += corr - 1.;
223     dsum2 += (corr-1.)*(corr-1.);
224 greg 2.3 }
225 greg 2.6 lastVar = thisVar;
226     thisVar = dsum2/(double)nn;
227 greg 2.10 #ifdef DEBUG
228 greg 2.4 fprintf(stderr, "Avg., RMS error: %.1f%% %.1f%%\n",
229     100.*dsum/(double)nn,
230 greg 2.6 100.*sqrt(thisVar));
231 greg 2.10 #endif
232 greg 2.6 } while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
233    
234 greg 2.8 nn = 0; /* compute sum for normalization */
235     while (nn < newnode->nrbf)
236     newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
237    
238 greg 2.10 insert_dsf(newnode);
239     return(newnode);
240 greg 2.3 }
241    
242 greg 2.1 /* Load a set of measurements corresponding to a particular incident angle */
243     static int
244 greg 2.10 load_pabopto_meas(const char *fname)
245 greg 2.1 {
246     FILE *fp = fopen(fname, "r");
247     int inp_is_DSF = -1;
248 greg 2.10 double new_phi, theta_out, phi_out, val;
249 greg 2.1 char buf[2048];
250     int n, c;
251    
252     if (fp == NULL) {
253     fputs(fname, stderr);
254     fputs(": cannot open\n", stderr);
255     return(0);
256     }
257 greg 2.5 memset(dsf_grid, 0, sizeof(dsf_grid));
258 greg 2.10 #ifdef DEBUG
259     fprintf(stderr, "Loading measurement file '%s'...\n", fname);
260     #endif
261 greg 2.1 /* read header information */
262     while ((c = getc(fp)) == '#' || c == EOF) {
263     if (fgets(buf, sizeof(buf), fp) == NULL) {
264     fputs(fname, stderr);
265     fputs(": unexpected EOF\n", stderr);
266     fclose(fp);
267     return(0);
268     }
269     if (!strcmp(buf, "format: theta phi DSF\n")) {
270     inp_is_DSF = 1;
271     continue;
272     }
273     if (!strcmp(buf, "format: theta phi BSDF\n")) {
274     inp_is_DSF = 0;
275     continue;
276     }
277     if (sscanf(buf, "intheta %lf", &theta_in_deg) == 1)
278     continue;
279 greg 2.10 if (sscanf(buf, "inphi %lf", &new_phi) == 1)
280 greg 2.1 continue;
281     if (sscanf(buf, "incident_angle %lf %lf",
282 greg 2.10 &theta_in_deg, &new_phi) == 2)
283 greg 2.1 continue;
284     }
285     if (inp_is_DSF < 0) {
286     fputs(fname, stderr);
287     fputs(": unknown format\n", stderr);
288     fclose(fp);
289     return(0);
290     }
291 greg 2.10 if (!input_orient) /* check input orientation */
292     input_orient = 1 - 2*(theta_in_deg > 90.);
293     else if (input_orient > 0 ^ theta_in_deg < 90.) {
294     fputs("Cannot handle input angles on both sides of surface\n",
295     stderr);
296     exit(1);
297     }
298     if (single_plane_incident > 0) /* check if still in plane */
299     single_plane_incident = (round(new_phi) == round(phi_in_deg));
300     else if (single_plane_incident < 0)
301     single_plane_incident = 1;
302     phi_in_deg = new_phi;
303     ungetc(c, fp); /* read actual data */
304 greg 2.1 while (fscanf(fp, "%lf %lf %lf\n", &theta_out, &phi_out, &val) == 3) {
305     FVECT ovec;
306     int pos[2];
307    
308 greg 2.10 if (!output_orient) /* check output orientation */
309     output_orient = 1 - 2*(theta_out > 90.);
310     else if (output_orient > 0 ^ theta_out < 90.) {
311     fputs("Cannot handle output angles on both sides of surface\n",
312     stderr);
313     exit(1);
314     }
315 greg 2.1 ovec[2] = sin(M_PI/180.*theta_out);
316     ovec[0] = cos(M_PI/180.*phi_out) * ovec[2];
317     ovec[1] = sin(M_PI/180.*phi_out) * ovec[2];
318     ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
319    
320 greg 2.5 if (!inp_is_DSF)
321     val *= ovec[2]; /* convert from BSDF to DSF */
322 greg 2.1
323     pos_from_vec(pos, ovec);
324    
325 greg 2.5 dsf_grid[pos[0]][pos[1]].vsum += val;
326     dsf_grid[pos[0]][pos[1]].nval++;
327 greg 2.1 }
328     n = 0;
329     while ((c = getc(fp)) != EOF)
330     n += !isspace(c);
331     if (n)
332     fprintf(stderr,
333     "%s: warning: %d unexpected characters past EOD\n",
334     fname, n);
335     fclose(fp);
336     return(1);
337     }
338    
339     /* Compute radii for non-empty bins */
340     /* (distance to furthest empty bin for which non-empty bin is the closest) */
341     static void
342     compute_radii(void)
343     {
344 greg 2.4 unsigned int fill_grid[GRIDRES][GRIDRES];
345     unsigned short fill_cnt[GRIDRES][GRIDRES];
346 greg 2.2 FVECT ovec0, ovec1;
347     double ang2, lastang2;
348     int r, i, j, jn, ii, jj, inear, jnear;
349    
350     r = GRIDRES/2; /* proceed in zig-zag */
351 greg 2.1 for (i = 0; i < GRIDRES; i++)
352     for (jn = 0; jn < GRIDRES; jn++) {
353     j = (i&1) ? jn : GRIDRES-1-jn;
354 greg 2.5 if (dsf_grid[i][j].nval) /* find empty grid pos. */
355 greg 2.1 continue;
356 greg 2.10 ovec_from_pos(ovec0, i, j);
357 greg 2.1 inear = jnear = -1; /* find nearest non-empty */
358 greg 2.2 lastang2 = M_PI*M_PI;
359 greg 2.1 for (ii = i-r; ii <= i+r; ii++) {
360     if (ii < 0) continue;
361     if (ii >= GRIDRES) break;
362     for (jj = j-r; jj <= j+r; jj++) {
363     if (jj < 0) continue;
364     if (jj >= GRIDRES) break;
365 greg 2.5 if (!dsf_grid[ii][jj].nval)
366 greg 2.1 continue;
367 greg 2.10 ovec_from_pos(ovec1, ii, jj);
368 greg 2.2 ang2 = 2. - 2.*DOT(ovec0,ovec1);
369     if (ang2 >= lastang2)
370 greg 2.1 continue;
371 greg 2.2 lastang2 = ang2;
372 greg 2.1 inear = ii; jnear = jj;
373     }
374     }
375 greg 2.2 if (inear < 0) {
376     fputs("Could not find non-empty neighbor!\n", stderr);
377     exit(1);
378     }
379     ang2 = sqrt(lastang2);
380     r = ANG2R(ang2); /* record if > previous */
381 greg 2.5 if (r > dsf_grid[inear][jnear].crad)
382     dsf_grid[inear][jnear].crad = r;
383 greg 2.2 /* next search radius */
384 greg 2.10 r = ang2*(2.*GRIDRES/M_PI) + 3;
385 greg 2.1 }
386 greg 2.4 /* blur radii over hemisphere */
387 greg 2.1 memset(fill_grid, 0, sizeof(fill_grid));
388 greg 2.4 memset(fill_cnt, 0, sizeof(fill_cnt));
389 greg 2.1 for (i = 0; i < GRIDRES; i++)
390     for (j = 0; j < GRIDRES; j++) {
391 greg 2.5 if (!dsf_grid[i][j].crad)
392 greg 2.4 continue; /* missing distance */
393 greg 2.5 r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
394 greg 2.1 for (ii = i-r; ii <= i+r; ii++) {
395     if (ii < 0) continue;
396     if (ii >= GRIDRES) break;
397     for (jj = j-r; jj <= j+r; jj++) {
398     if (jj < 0) continue;
399     if (jj >= GRIDRES) break;
400 greg 2.4 if ((ii-i)*(ii-i) + (jj-j)*(jj-j) > r*r)
401 greg 2.1 continue;
402 greg 2.5 fill_grid[ii][jj] += dsf_grid[i][j].crad;
403 greg 2.4 fill_cnt[ii][jj]++;
404 greg 2.1 }
405     }
406     }
407 greg 2.6 /* copy back blurred radii */
408 greg 2.1 for (i = 0; i < GRIDRES; i++)
409     for (j = 0; j < GRIDRES; j++)
410 greg 2.4 if (fill_cnt[i][j])
411 greg 2.5 dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
412 greg 2.1 }
413    
414 greg 2.6 /* Cull points for more uniform distribution, leave all nval 0 or 1 */
415 greg 2.1 static void
416     cull_values(void)
417     {
418 greg 2.2 FVECT ovec0, ovec1;
419     double maxang, maxang2;
420     int i, j, ii, jj, r;
421 greg 2.1 /* simple greedy algorithm */
422     for (i = 0; i < GRIDRES; i++)
423     for (j = 0; j < GRIDRES; j++) {
424 greg 2.5 if (!dsf_grid[i][j].nval)
425 greg 2.1 continue;
426 greg 2.5 if (!dsf_grid[i][j].crad)
427 greg 2.2 continue; /* shouldn't happen */
428 greg 2.10 ovec_from_pos(ovec0, i, j);
429 greg 2.5 maxang = 2.*R2ANG(dsf_grid[i][j].crad);
430 greg 2.2 if (maxang > ovec0[2]) /* clamp near horizon */
431     maxang = ovec0[2];
432     r = maxang*(2.*GRIDRES/M_PI) + 1;
433     maxang2 = maxang*maxang;
434 greg 2.1 for (ii = i-r; ii <= i+r; ii++) {
435     if (ii < 0) continue;
436     if (ii >= GRIDRES) break;
437     for (jj = j-r; jj <= j+r; jj++) {
438     if (jj < 0) continue;
439     if (jj >= GRIDRES) break;
440 greg 2.5 if (!dsf_grid[ii][jj].nval)
441 greg 2.1 continue;
442 greg 2.2 if ((ii == i) & (jj == j))
443     continue; /* don't get self-absorbed */
444 greg 2.10 ovec_from_pos(ovec1, ii, jj);
445 greg 2.2 if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
446 greg 2.1 continue;
447 greg 2.2 /* absorb sum */
448 greg 2.5 dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
449     dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
450 greg 2.2 /* keep value, though */
451 greg 2.6 dsf_grid[ii][jj].vsum /= (float)dsf_grid[ii][jj].nval;
452 greg 2.5 dsf_grid[ii][jj].nval = 0;
453 greg 2.1 }
454     }
455     }
456 greg 2.6 /* final averaging pass */
457     for (i = 0; i < GRIDRES; i++)
458     for (j = 0; j < GRIDRES; j++)
459     if (dsf_grid[i][j].nval > 1) {
460     dsf_grid[i][j].vsum /= (float)dsf_grid[i][j].nval;
461     dsf_grid[i][j].nval = 1;
462     }
463 greg 2.1 }
464    
465 greg 2.8 /* Compute (and allocate) migration price matrix for optimization */
466     static float *
467 greg 2.10 price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
468 greg 2.8 {
469     float *pmtx = (float *)malloc(sizeof(float) *
470     from_rbf->nrbf * to_rbf->nrbf);
471     FVECT *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
472     int i, j;
473    
474     if ((pmtx == NULL) | (vto == NULL)) {
475     fputs("Out of memory in migration_costs()\n", stderr);
476     exit(1);
477     }
478     for (j = to_rbf->nrbf; j--; ) /* save repetitive ops. */
479 greg 2.10 ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
480 greg 2.8
481     for (i = from_rbf->nrbf; i--; ) {
482     const double from_ang = R2ANG(from_rbf->rbfa[i].crad);
483     FVECT vfrom;
484 greg 2.10 ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
485 greg 2.8 for (j = to_rbf->nrbf; j--; )
486     pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
487     fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
488     }
489     free(vto);
490     return(pmtx);
491     }
492    
493     /* Comparison routine needed for sorting price row */
494     static const float *price_arr;
495     static int
496     msrt_cmp(const void *p1, const void *p2)
497     {
498     float c1 = price_arr[*(const int *)p1];
499     float c2 = price_arr[*(const int *)p2];
500    
501     if (c1 > c2) return(1);
502     if (c1 < c2) return(-1);
503     return(0);
504     }
505    
506     /* Compute minimum (optimistic) cost for moving the given source material */
507     static double
508     min_cost(double amt2move, const double *avail, const float *price, int n)
509     {
510     static int *price_sort = NULL;
511     static int n_alloc = 0;
512     double total_cost = 0;
513     int i;
514    
515     if (amt2move <= FTINY) /* pre-emptive check */
516     return(0.);
517     if (n > n_alloc) { /* (re)allocate sort array */
518     if (n_alloc) free(price_sort);
519     price_sort = (int *)malloc(sizeof(int)*n);
520     if (price_sort == NULL) {
521     fputs("Out of memory in min_cost()\n", stderr);
522     exit(1);
523     }
524     n_alloc = n;
525     }
526     for (i = n; i--; )
527     price_sort[i] = i;
528     price_arr = price;
529     qsort(price_sort, n, sizeof(int), &msrt_cmp);
530     /* move cheapest first */
531     for (i = 0; i < n && amt2move > FTINY; i++) {
532     int d = price_sort[i];
533     double amt = (amt2move < avail[d]) ? amt2move : avail[d];
534    
535     total_cost += amt * price[d];
536     amt2move -= amt;
537     }
538     return(total_cost);
539     }
540    
541     /* Take a step in migration by choosing optimal bucket to transfer */
542     static double
543     migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
544     {
545     static double *src_cost = NULL;
546     int n_alloc = 0;
547 greg 2.9 const double maxamt = 2./(mtx_nrows(mig)*mtx_ncols(mig));
548 greg 2.8 double amt = 0;
549     struct {
550     int s, d; /* source and destination */
551     double price; /* price estimate per amount moved */
552     double amt; /* amount we can move */
553     } cur, best;
554     int i;
555    
556     if (mtx_nrows(mig) > n_alloc) { /* allocate cost array */
557     if (n_alloc)
558     free(src_cost);
559     src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
560     if (src_cost == NULL) {
561     fputs("Out of memory in migration_step()\n", stderr);
562     exit(1);
563     }
564     n_alloc = mtx_nrows(mig);
565     }
566     for (i = mtx_nrows(mig); i--; ) /* starting costs for diff. */
567     src_cost[i] = min_cost(src_rem[i], dst_rem,
568     pmtx+i*mtx_ncols(mig), mtx_ncols(mig));
569    
570     /* find best source & dest. */
571     best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
572     for (cur.s = mtx_nrows(mig); cur.s--; ) {
573     const float *price = pmtx + cur.s*mtx_ncols(mig);
574     double cost_others = 0;
575     if (src_rem[cur.s] <= FTINY)
576     continue;
577     cur.d = -1; /* examine cheapest dest. */
578     for (i = mtx_ncols(mig); i--; )
579     if (dst_rem[i] > FTINY &&
580     (cur.d < 0 || price[i] < price[cur.d]))
581     cur.d = i;
582     if (cur.d < 0)
583     return(.0);
584     if ((cur.price = price[cur.d]) >= best.price)
585     continue; /* no point checking further */
586     cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
587     src_rem[cur.s] : dst_rem[cur.d];
588     if (cur.amt > maxamt) cur.amt = maxamt;
589     dst_rem[cur.d] -= cur.amt; /* add up differential costs */
590     for (i = mtx_nrows(mig); i--; ) {
591     if (i == cur.s) continue;
592     cost_others += min_cost(src_rem[i], dst_rem, price, mtx_ncols(mig))
593     - src_cost[i];
594     }
595     dst_rem[cur.d] += cur.amt; /* undo trial move */
596     cur.price += cost_others/cur.amt; /* adjust effective price */
597     if (cur.price < best.price) /* are we better than best? */
598     best = cur;
599     }
600     if ((best.s < 0) | (best.d < 0))
601     return(.0);
602     /* make the actual move */
603     mig->mtx[mtx_ndx(mig,best.s,best.d)] += best.amt;
604     src_rem[best.s] -= best.amt;
605     dst_rem[best.d] -= best.amt;
606     return(best.amt);
607     }
608    
609     /* Compute (and insert) migration along directed edge */
610     static MIGRATION *
611 greg 2.10 make_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
612 greg 2.8 {
613     const double end_thresh = 0.02/(from_rbf->nrbf*to_rbf->nrbf);
614     float *pmtx = price_routes(from_rbf, to_rbf);
615     MIGRATION *newmig = (MIGRATION *)malloc(sizeof(MIGRATION) +
616     sizeof(float) *
617     (from_rbf->nrbf*to_rbf->nrbf - 1));
618     double *src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
619     double *dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
620     double total_rem = 1.;
621     int i;
622    
623     if ((newmig == NULL) | (src_rem == NULL) | (dst_rem == NULL)) {
624     fputs("Out of memory in make_migration()\n", stderr);
625     exit(1);
626     }
627 greg 2.10 #ifdef DEBUG
628     {
629     double theta, phi;
630     theta = acos(from_rbf->invec[2])*(180./M_PI);
631     phi = atan2(from_rbf->invec[1],from_rbf->invec[0])*(180./M_PI);
632     fprintf(stderr, "Building path from (theta,phi) (%d,%d) to ",
633     round(theta), round(phi));
634     theta = acos(to_rbf->invec[2])*(180./M_PI);
635     phi = atan2(to_rbf->invec[1],to_rbf->invec[0])*(180./M_PI);
636     fprintf(stderr, "(%d,%d)\n", round(theta), round(phi));
637     }
638     #endif
639 greg 2.8 newmig->next = NULL;
640     newmig->rbfv[0] = from_rbf;
641     newmig->rbfv[1] = to_rbf;
642     newmig->enxt[0] = newmig->enxt[1] = NULL;
643     memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
644     /* starting quantities */
645     for (i = from_rbf->nrbf; i--; )
646     src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
647     for (i = to_rbf->nrbf; i--; )
648     dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
649     /* move a bit at a time */
650     while (total_rem > end_thresh)
651     total_rem -= migration_step(newmig, src_rem, dst_rem, pmtx);
652    
653     free(pmtx); /* free working arrays */
654     free(src_rem);
655     free(dst_rem);
656     for (i = from_rbf->nrbf; i--; ) { /* normalize final matrix */
657     float nf = rbf_volume(&from_rbf->rbfa[i]);
658     int j;
659     if (nf <= FTINY) continue;
660     nf = from_rbf->vtotal / nf;
661     for (j = to_rbf->nrbf; j--; )
662     newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
663     }
664     /* insert in edge lists */
665     newmig->enxt[0] = from_rbf->ejl;
666     from_rbf->ejl = newmig;
667     newmig->enxt[1] = to_rbf->ejl;
668     to_rbf->ejl = newmig;
669     newmig->next = mig_list;
670     return(mig_list = newmig);
671     }
672    
673 greg 2.10 /* Get triangle surface orientation (unnormalized) */
674     static void
675     tri_orient(FVECT vres, const FVECT v1, const FVECT v2, const FVECT v3)
676     {
677     FVECT v2minus1, v3minus2;
678    
679     VSUB(v2minus1, v2, v1);
680     VSUB(v3minus2, v3, v2);
681     VCROSS(vres, v2minus1, v3minus2);
682     }
683    
684     /* Determine if vertex order is reversed (inward normal) */
685     static int
686     is_rev_tri(const FVECT v1, const FVECT v2, const FVECT v3)
687     {
688     FVECT tor;
689    
690     tri_orient(tor, v1, v2, v3);
691    
692     return(DOT(tor, v2) < 0.);
693     }
694    
695     /* Find vertices completing triangles on either side of the given edge */
696     static int
697     get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
698     {
699     const MIGRATION *ej, *ej2;
700     RBFNODE *tv;
701    
702     rbfv[0] = rbfv[1] = NULL;
703     for (ej = mig->rbfv[0]->ejl; ej != NULL;
704     ej = nextedge(mig->rbfv[0],ej)) {
705     if (ej == mig)
706     continue;
707     tv = opp_rbf(mig->rbfv[0],ej);
708     for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
709     if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
710     rbfv[is_rev_tri(mig->rbfv[0]->invec,
711     mig->rbfv[1]->invec,
712     tv->invec)] = tv;
713     break;
714     }
715     }
716     return((rbfv[0] != NULL) + (rbfv[1] != NULL));
717     }
718    
719     /* Find context hull vertex to complete triangle (oriented call) */
720     static RBFNODE *
721     find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
722 greg 2.8 {
723 greg 2.10 FVECT vmid, vor;
724     RBFNODE *rbf, *rbfbest = NULL;
725     double dprod2, bestdprod2 = 0.5;
726    
727     VADD(vmid, rbf0->invec, rbf1->invec);
728     if (normalize(vmid) == 0)
729     return(NULL);
730     /* XXX exhaustive search */
731     for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
732     if ((rbf == rbf0) | (rbf == rbf1))
733     continue;
734     tri_orient(vor, rbf0->invec, rbf1->invec, rbf->invec);
735     dprod2 = DOT(vor, vmid);
736     if (dprod2 <= FTINY)
737     continue; /* wrong orientation */
738     dprod2 *= dprod2 / DOT(vor,vor);
739     if (dprod2 > bestdprod2) { /* more convex than prev? */
740     rbfbest = rbf;
741     bestdprod2 = dprod2;
742     }
743     }
744     return(rbf);
745     }
746    
747     /* Create new migration edge and grow mesh recursively around it */
748     static void
749     mesh_from_edge(RBFNODE *rbf0, RBFNODE *rbf1)
750     {
751     MIGRATION *newej;
752     RBFNODE *tvert[2];
753    
754     if (rbf0 < rbf1) /* avoid migration loops */
755     newej = make_migration(rbf0, rbf1);
756     else
757     newej = make_migration(rbf1, rbf0);
758     /* triangle on either side? */
759     get_triangles(tvert, newej);
760     if (tvert[0] == NULL) { /* recurse on new right edge */
761     tvert[0] = find_chull_vert(newej->rbfv[0], newej->rbfv[1]);
762     if (tvert[0] != NULL) {
763     mesh_from_edge(rbf0, tvert[0]);
764     mesh_from_edge(rbf1, tvert[0]);
765     }
766     }
767     if (tvert[1] == NULL) { /* recurse on new left edge */
768     tvert[1] = find_chull_vert(newej->rbfv[1], newej->rbfv[0]);
769     if (tvert[1] != NULL) {
770     mesh_from_edge(rbf0, tvert[1]);
771     mesh_from_edge(rbf1, tvert[1]);
772     }
773     }
774     }
775 greg 2.8
776 greg 2.10 /* Draw edge list into mig_grid array */
777     static void
778     draw_edges()
779     {
780     int nnull = 0, ntot = 0;
781     MIGRATION *ej;
782     int p0[2], p1[2];
783    
784     /* memset(mig_grid, 0, sizeof(mig_grid)); */
785     for (ej = mig_list; ej != NULL; ej = ej->next) {
786     ++ntot;
787     pos_from_vec(p0, ej->rbfv[0]->invec);
788     pos_from_vec(p1, ej->rbfv[1]->invec);
789     if ((p0[0] == p1[0]) & (p0[1] == p1[1])) {
790     ++nnull;
791     mig_grid[p0[0]][p0[1]] = ej;
792     continue;
793     }
794     if (abs(p1[0]-p0[0]) > abs(p1[1]-p0[1])) {
795     const int xstep = 2*(p1[0] > p0[0]) - 1;
796     const double ystep = (double)((p1[1]-p0[1])*xstep) /
797     (double)(p1[0]-p0[0]);
798     int x;
799     double y;
800     for (x = p0[0], y = p0[1]+.5; x != p1[0];
801     x += xstep, y += ystep)
802     mig_grid[x][(int)y] = ej;
803     mig_grid[x][(int)y] = ej;
804     } else {
805     const int ystep = 2*(p1[1] > p0[1]) - 1;
806     const double xstep = (double)((p1[0]-p0[0])*ystep) /
807     (double)(p1[1]-p0[1]);
808     int y;
809     double x;
810     for (y = p0[1], x = p0[0]+.5; y != p1[1];
811     y += ystep, x += xstep)
812     mig_grid[(int)x][y] = ej;
813     mig_grid[(int)x][y] = ej;
814     }
815     }
816     if (nnull)
817     fprintf(stderr, "Warning: %d of %d edges are null\n",
818     nnull, ntot);
819     }
820    
821     /* Build our triangle mesh from recorded RBFs */
822     static void
823     build_mesh()
824     {
825     double best2 = M_PI*M_PI;
826     RBFNODE *rbf, *rbf_near = NULL;
827     /* check if isotropic */
828     if (single_plane_incident) {
829     for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
830     if (rbf->next != NULL)
831     make_migration(rbf, rbf->next);
832     return;
833     }
834     /* find RBF nearest to head */
835     if (dsf_list == NULL)
836     return;
837     for (rbf = dsf_list->next; rbf != NULL; rbf = rbf->next) {
838     double dist2 = 2. - 2.*DOT(dsf_list->invec,rbf->invec);
839     if (dist2 < best2) {
840     rbf_near = rbf;
841     best2 = dist2;
842     }
843     }
844     if (rbf_near == NULL) {
845     fputs("Cannot find nearest point for first edge\n", stderr);
846 greg 2.8 exit(1);
847     }
848 greg 2.10 /* build mesh from this edge */
849     mesh_from_edge(dsf_list, rbf_near);
850     /* draw edge list into grid */
851     draw_edges();
852     }
853    
854     /* Identify enclosing triangle for this position (flood fill raster check) */
855     static int
856     identify_tri(MIGRATION *miga[3], unsigned char vmap[GRIDRES][(GRIDRES+7)/8],
857     int px, int py)
858     {
859     const int btest = 1<<(py&07);
860    
861     if (vmap[px][py>>3] & btest) /* already visited here? */
862     return(1);
863     /* else mark it */
864     vmap[px][py>>3] |= btest;
865    
866     if (mig_grid[px][py] != NULL) { /* are we on an edge? */
867     int i;
868     for (i = 0; i < 3; i++) {
869     if (miga[i] == mig_grid[px][py])
870     return(1);
871     if (miga[i] != NULL)
872     continue;
873     while (i > 0 && miga[i-1] > mig_grid[px][py]) {
874     miga[i] = miga[i-1];
875     --i; /* order vertices by pointer */
876     }
877     miga[i] = mig_grid[px][py];
878     return(1);
879     }
880     return(0); /* outside triangle! */
881     }
882     /* check neighbors (flood) */
883     if (px > 0 && !identify_tri(miga, vmap, px-1, py))
884     return(0);
885     if (px < GRIDRES-1 && !identify_tri(miga, vmap, px+1, py))
886     return(0);
887     if (py > 0 && !identify_tri(miga, vmap, px, py-1))
888     return(0);
889     if (py < GRIDRES-1 && !identify_tri(miga, vmap, px, py+1))
890     return(0);
891     return(1); /* this neighborhood done */
892     }
893    
894     /* Find edge(s) for interpolating the given incident vector */
895     static int
896     get_interp(MIGRATION *miga[3], const FVECT invec)
897     {
898     miga[0] = miga[1] = miga[2] = NULL;
899     if (single_plane_incident) { /* isotropic BSDF? */
900     RBFNODE *rbf; /* find edge we're on */
901     for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
902     if (input_orient*rbf->invec[2] < input_orient*invec[2])
903     break;
904     if (rbf->next != NULL &&
905     input_orient*rbf->next->invec[2] <
906     input_orient*invec[2]) {
907     for (miga[0] = rbf->ejl; miga[0] != NULL;
908     miga[0] = nextedge(rbf,miga[0]))
909     if (opp_rbf(rbf,miga[0]) == rbf->next)
910     return(1);
911     break;
912     }
913     }
914     return(0); /* outside range! */
915     }
916     { /* else use triagnle mesh */
917     unsigned char floodmap[GRIDRES][(GRIDRES+7)/8];
918     int pstart[2];
919    
920     pos_from_vec(pstart, invec);
921     memset(floodmap, 0, sizeof(floodmap));
922     /* call flooding function */
923     if (!identify_tri(miga, floodmap, pstart[0], pstart[1]))
924     return(0); /* outside mesh */
925     if ((miga[0] == NULL) | (miga[2] == NULL))
926     return(0); /* should never happen */
927     if (miga[1] == NULL)
928     return(1); /* on edge */
929     return(3); /* else in triangle */
930     }
931     }
932    
933     /* Advect and allocate new RBF along edge */
934     static RBFNODE *
935     e_advect_rbf(const MIGRATION *mig, const FVECT invec)
936     {
937     RBFNODE *rbf;
938     int n, i, j;
939     double t, full_dist;
940     /* get relative position */
941     t = acos(DOT(invec, mig->rbfv[0]->invec));
942     if (t < M_PI/GRIDRES) { /* near first DSF */
943     n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
944     rbf = (RBFNODE *)malloc(n);
945     if (rbf == NULL)
946     goto memerr;
947     memcpy(rbf, mig->rbfv[0], n); /* just duplicate */
948     return(rbf);
949     }
950     full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
951     if (t > full_dist-M_PI/GRIDRES) { /* near second DSF */
952     n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
953     rbf = (RBFNODE *)malloc(n);
954     if (rbf == NULL)
955     goto memerr;
956     memcpy(rbf, mig->rbfv[1], n); /* just duplicate */
957     return(rbf);
958     }
959     t /= full_dist;
960     n = 0; /* count migrating particles */
961     for (i = 0; i < mtx_nrows(mig); i++)
962     for (j = 0; j < mtx_ncols(mig); j++)
963     n += (mig->mtx[mtx_ndx(mig,i,j)] > FTINY);
964    
965     rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
966     if (rbf == NULL)
967     goto memerr;
968     rbf->next = NULL; rbf->ejl = NULL;
969     VCOPY(rbf->invec, invec);
970     rbf->nrbf = n;
971     rbf->vtotal = 1.-t + t*mig->rbfv[1]->vtotal/mig->rbfv[0]->vtotal;
972     n = 0; /* advect RBF lobes */
973     for (i = 0; i < mtx_nrows(mig); i++) {
974     const RBFVAL *rbf0i = &mig->rbfv[0]->rbfa[i];
975     const float peak0 = rbf0i->peak;
976     const double rad0 = R2ANG(rbf0i->crad);
977     FVECT v0;
978     float mv;
979     ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
980     for (j = 0; j < mtx_ncols(mig); j++)
981     if ((mv = mig->mtx[mtx_ndx(mig,i,j)]) > FTINY) {
982     const RBFVAL *rbf1j = &mig->rbfv[1]->rbfa[j];
983     double rad1 = R2ANG(rbf1j->crad);
984     FVECT v;
985     int pos[2];
986     rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal;
987     rbf->rbfa[n].crad = ANG2R(sqrt(rad0*rad0*(1.-t) +
988     rad1*rad1*t));
989     ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
990     geodesic(v, v0, v, t, GEOD_REL);
991     pos_from_vec(pos, v);
992     rbf->rbfa[n].gx = pos[0];
993     rbf->rbfa[n].gy = pos[1];
994     ++n;
995     }
996     }
997     rbf->vtotal *= mig->rbfv[0]->vtotal; /* turn ratio into actual */
998     return(rbf);
999     memerr:
1000     fputs("Out of memory in e_advect_rbf()\n", stderr);
1001     exit(1);
1002     return(NULL); /* pro forma return */
1003     }
1004    
1005     /* Partially advect between recorded incident angles and allocate new RBF */
1006     static RBFNODE *
1007     advect_rbf(const FVECT invec)
1008     {
1009     MIGRATION *miga[3];
1010     RBFNODE *rbf;
1011     int n, i, j;
1012     double s, t;
1013    
1014     if (!get_interp(miga, invec)) /* can't interpolate? */
1015     return(NULL);
1016     if (miga[1] == NULL) /* along edge? */
1017     return(e_advect_rbf(miga[0], invec));
1018     /* figure out position */
1019    
1020     rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
1021 greg 2.8 if (rbf == NULL) {
1022     fputs("Out of memory in advect_rbf()\n", stderr);
1023     exit(1);
1024     }
1025     rbf->next = NULL; rbf->ejl = NULL;
1026     VCOPY(rbf->invec, invec);
1027     rbf->vtotal = 0;
1028 greg 2.10 rbf->nrbf = n;
1029 greg 2.8
1030 greg 2.10 return(rbf);
1031 greg 2.8 }
1032 greg 2.1
1033     #if 1
1034     /* Test main produces a Radiance model from the given input file */
1035     int
1036     main(int argc, char *argv[])
1037     {
1038     char buf[128];
1039     FILE *pfp;
1040     double bsdf;
1041     FVECT dir;
1042     int i, j, n;
1043    
1044     if (argc != 2) {
1045     fprintf(stderr, "Usage: %s input.dat > output.rad\n", argv[0]);
1046     return(1);
1047     }
1048 greg 2.10 if (!load_pabopto_meas(argv[1]))
1049 greg 2.1 return(1);
1050    
1051     compute_radii();
1052     cull_values();
1053 greg 2.3 make_rbfrep();
1054     /* produce spheres at meas. */
1055     puts("void plastic yellow\n0\n0\n5 .6 .4 .01 .04 .08\n");
1056 greg 2.1 puts("void plastic pink\n0\n0\n5 .5 .05 .9 .04 .08\n");
1057     n = 0;
1058     for (i = 0; i < GRIDRES; i++)
1059     for (j = 0; j < GRIDRES; j++)
1060 greg 2.5 if (dsf_grid[i][j].vsum > .0f) {
1061 greg 2.10 ovec_from_pos(dir, i, j);
1062 greg 2.5 bsdf = dsf_grid[i][j].vsum / dir[2];
1063     if (dsf_grid[i][j].nval) {
1064 greg 2.3 printf("pink cone c%04d\n0\n0\n8\n", ++n);
1065     printf("\t%.6g %.6g %.6g\n",
1066 greg 2.1 dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
1067 greg 2.3 printf("\t%.6g %.6g %.6g\n",
1068 greg 2.1 dir[0]*(bsdf+.005), dir[1]*(bsdf+.005),
1069     dir[2]*(bsdf+.005));
1070 greg 2.3 puts("\t.003\t0\n");
1071     } else {
1072 greg 2.10 ovec_from_pos(dir, i, j);
1073 greg 2.3 printf("yellow sphere s%04d\n0\n0\n", ++n);
1074     printf("4 %.6g %.6g %.6g .0015\n\n",
1075     dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
1076     }
1077 greg 2.1 }
1078     /* output continuous surface */
1079     puts("void trans tgreen\n0\n0\n7 .7 1 .7 .04 .04 .9 .9\n");
1080     fflush(stdout);
1081 greg 2.5 sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES-1, GRIDRES-1);
1082 greg 2.1 pfp = popen(buf, "w");
1083     if (pfp == NULL) {
1084     fputs(buf, stderr);
1085     fputs(": cannot start command\n", stderr);
1086     return(1);
1087     }
1088     for (i = 0; i < GRIDRES; i++)
1089     for (j = 0; j < GRIDRES; j++) {
1090 greg 2.10 ovec_from_pos(dir, i, j);
1091 greg 2.5 bsdf = eval_rbfrep(dsf_list, dir) / dir[2];
1092 greg 2.1 fprintf(pfp, "%.8e %.8e %.8e\n",
1093     dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
1094     }
1095     return(pclose(pfp)==0 ? 0 : 1);
1096     }
1097     #endif