ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rhoptimize.c
Revision: 3.2
Committed: Thu Nov 5 14:43:14 1998 UTC (25 years, 4 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.1: +63 -28 lines
Log Message:
build queue then read
bug fixes in neighbor algorithm

File Contents

# User Rev Content
1 gwlarson 3.1 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2    
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ SGI";
5     #endif
6    
7     /*
8     * Optimize holodeck for quick access.
9     *
10     * 11/4/98 Greg Ward Larson
11     */
12    
13     #include "holo.h"
14    
15     #ifndef BKBSIZE
16 gwlarson 3.2 #define BKBSIZE 256 /* beam clump size (kilobytes) */
17 gwlarson 3.1 #endif
18    
19 gwlarson 3.2 #define flgop(p,i,op) ((p)[(i)>>5] op (1L<<((i)&0x1f)))
20     #define isset(p,i) flgop(p,i,&)
21     #define setfl(p,i) flgop(p,i,|=)
22     #define clrfl(p,i) flgop(p,i,&=~)
23 gwlarson 3.1
24     char *progname;
25    
26     extern long rhinitcopy();
27    
28    
29     main(argc, argv)
30     int argc;
31     char *argv[];
32     {
33     char nambuf[128];
34     char *inpname, *outname;
35     int hdfd[2];
36     long nextipos, lastopos, thisopos;
37    
38     progname = argv[0];
39     if (argc < 2 | argc > 3) {
40     fprintf(stderr, "Usage: %s input.hdk [output.hdk]\n", progname);
41     exit(1);
42     }
43     inpname = argv[1];
44     if (argc == 3) /* use given output file */
45     outname = argv[2];
46     else { /* else create temporary file */
47     strcpy(nambuf, inpname);
48     if ((outname = strrchr(nambuf, '/')) != NULL)
49     outname++;
50     else
51     outname = nambuf;
52     sprintf(outname, "rho%d.hdk", getpid());
53     outname = nambuf;
54     }
55     /* copy holodeck file header */
56     nextipos = rhinitcopy(hdfd, inpname, outname);
57     lastopos = 0L; /* copy sections one by one */
58     while (nextipos != 0L) {
59     /* set input position; get next */
60     lseek(hdfd[0], nextipos, 0);
61     read(hdfd[0], (char *)&nextipos, sizeof(nextipos));
62     /* get output position; set last */
63     thisopos = lseek(hdfd[1], 0L, 2);
64     if (lastopos > 0L) {
65     lseek(hdfd[1], lastopos, 0);
66     write(hdfd[1], (char *)&thisopos, sizeof(thisopos));
67     lseek(hdfd[1], 0L, 2);
68     }
69     lastopos = thisopos;
70     thisopos = 0L; /* write place holder */
71     write(hdfd[1], (char *)&thisopos, sizeof(thisopos));
72     /* copy holodeck section */
73     copysect(hdfd[0], hdfd[1]);
74     }
75     /* clean up */
76     close(hdfd[0]);
77     close(hdfd[1]);
78     if (argc == 2 && rename(outname, inpname) < 0) {
79     sprintf(errmsg, "cannot rename \"%s\" to \"%s\"",
80     outname, inpname);
81     error(SYSTEM, errmsg);
82     }
83     exit(0);
84     }
85    
86    
87     long
88     rhinitcopy(hfd, infn, outfn) /* open files and copy header */
89     int hfd[2]; /* returned file descriptors */
90     char *infn, *outfn;
91     {
92     FILE *infp, *outfp;
93     long ifpos;
94     /* open files for i/o */
95     if ((infp = fopen(infn, "r")) == NULL) {
96     sprintf(errmsg, "cannot open \"%s\" for reading", infn);
97     error(SYSTEM, errmsg);
98     }
99     if ((outfp = fopen(outfn, "w+")) == NULL) {
100     sprintf(errmsg, "cannot open \"%s\" for writing", outfn);
101     error(SYSTEM, errmsg);
102     }
103     /* copy and verify header */
104     if (checkheader(infp, HOLOFMT, outfp) < 0 ||
105     getw(infp) != HOLOMAGIC)
106     error(USER, "input not in holodeck format");
107     fputformat(HOLOFMT, outfp);
108     fputc('\n', outfp);
109     putw(HOLOMAGIC, outfp);
110     /* get descriptors and free stdio */
111     if ((hfd[0] = dup(fileno(infp))) < 0 ||
112     (hfd[1] = dup(fileno(outfp))) < 0)
113     error(SYSTEM, "dup call failed in rhinitcopy");
114     ifpos = ftell(infp);
115     fclose(infp);
116     if (fclose(outfp) == EOF)
117     error(SYSTEM, "file flushing error in rhinitcopy");
118     /* we flush everything manually */
119     hdcachesize = 0;
120     /* return input position */
121     return(ifpos);
122     }
123    
124    
125     gcshifti(gc, ia, di, hp) /* shift cell row or column */
126     register GCOORD *gc;
127     int ia, di;
128     register HOLO *hp;
129     {
130     int nw;
131    
132     if (di > 0) {
133     if (++gc->i[ia] >= hp->grid[((gc->w>>1)+1+ia)%3]) {
134     nw = ((gc->w&~1) + (ia<<1) + 3) % 6;
135     gc->i[ia] = gc->i[1-ia];
136     gc->i[1-ia] = gc->w&1 ? hp->grid[((nw>>1)+2-ia)%3]-1 : 0;
137     gc->w = nw;
138     }
139     } else if (di < 0) {
140     if (--gc->i[ia] < 0) {
141     nw = ((gc->w&~1) + (ia<<1) + 2) % 6;
142     gc->i[ia] = gc->i[1-ia];
143     gc->i[1-ia] = gc->w&1 ? hp->grid[((nw>>1)+2-ia)%3]-1 : 0;
144     gc->w = nw;
145     }
146     }
147     }
148    
149    
150     mkneighgrid(ng, hp, gc) /* compute neighborhood for grid cell */
151     GCOORD ng[3*3];
152     HOLO *hp;
153     GCOORD *gc;
154     {
155     GCOORD gci0;
156 gwlarson 3.2 register int i, j;
157 gwlarson 3.1
158     for (i = 3; i--; ) {
159     copystruct(&gci0, gc);
160     gcshifti(&gci0, 0, i-1, hp);
161     for (j = 3; j--; ) {
162     copystruct(ng+(3*i+j), &gci0);
163 gwlarson 3.2 gcshifti(ng+(3*i+j), gci0.w==gc->w, j-1, hp);
164 gwlarson 3.1 }
165     }
166     }
167    
168    
169     int bneighlist[9*9-1];
170     int bneighrem;
171    
172     #define nextneigh() (bneighrem<=0 ? 0 : bneighlist[--bneighrem])
173    
174     int
175     firstneigh(hp, b) /* initialize neighbor list and return first */
176     HOLO *hp;
177     int b;
178     {
179     GCOORD wg0[9], wg1[9], bgc[2];
180     int i, j;
181    
182     hdbcoord(bgc, hp, b);
183     mkneighgrid(wg0, hp, bgc);
184     mkneighgrid(wg1, hp, bgc+1);
185     bneighrem = 0;
186     for (i = 9; i--; )
187     for (j = 9; j--; ) {
188     if (i == 4 & j == 4)
189     continue;
190 gwlarson 3.2 if (wg0[i].w == wg1[j].w)
191     continue;
192 gwlarson 3.1 copystruct(bgc, wg0+i);
193     copystruct(bgc+1, wg1+j);
194     bneighlist[bneighrem++] = hdbindex(hp, bgc);
195 gwlarson 3.2 #ifdef DEBUG
196     if (bneighlist[bneighrem-1] <= 0)
197     error(CONSISTENCY, "bad beam in firstneigh");
198     #endif
199 gwlarson 3.1 }
200     return(nextneigh());
201     }
202    
203    
204 gwlarson 3.2 BEAMI *beamdir;
205    
206     int
207     bpcmp(b1p, b2p) /* compare beam positions on disk */
208     int *b1p, *b2p;
209     {
210     register long pdif = beamdir[*b1p].fo - beamdir[*b2p].fo;
211    
212     if (pdif > 0) return(1);
213     if (pdif < 0) return(-1);
214     return(0);
215     }
216    
217    
218 gwlarson 3.1 copysect(ifd, ofd) /* copy holodeck section from ifd to ofd */
219     int ifd, ofd;
220     {
221     static short primes[] = {9431,6803,4177,2659,1609,887,587,251,47,1};
222 gwlarson 3.2 register HOLO *hinp;
223     HOLO *hout;
224 gwlarson 3.1 register BEAM *bp;
225 gwlarson 3.2 unsigned int4 *bflags;
226 gwlarson 3.1 int *bqueue;
227     int bqlen;
228     int4 bqtotal;
229 gwlarson 3.2 int bc, bci, bqc, myprime;
230 gwlarson 3.1 register int i;
231     /* load input section directory */
232     hinp = hdinit(ifd, NULL);
233     /* create output section directory */
234     hout = hdinit(ofd, (HDGRID *)hinp);
235     /* allocate beam queue */
236 gwlarson 3.2 bqueue = (int *)malloc(nbeams(hinp)*sizeof(int));
237     bflags = (int4 *)calloc((nbeams(hinp)>>3)+1, sizeof(int4));
238     if (bqueue == NULL | bflags == NULL)
239 gwlarson 3.1 error(SYSTEM, "out of memory in copysect");
240 gwlarson 3.2 /* mark empty beams as done */
241     for (i = nbeams(hinp); i-- > 0; )
242     if (!hinp->bi[i].nrd)
243     setfl(bflags, i);
244 gwlarson 3.1 /* pick a good prime step size */
245     for (i = 0; primes[i]<<5 >= nbeams(hinp); i++)
246     ;
247     while ((myprime = primes[i++]) > 1)
248     if (nbeams(hinp) % myprime)
249     break;
250     /* add each input beam and neighbors */
251     for (bc = bci = nbeams(hinp); bc > 0; bc--,
252     bci += bci>myprime ? -myprime : nbeams(hinp)-myprime) {
253 gwlarson 3.2 if (isset(bflags, bci))
254 gwlarson 3.1 continue;
255     bqueue[0] = bci; /* initialize queue */
256     bqlen = 1;
257     bqtotal = bnrays(hinp, bci);
258 gwlarson 3.2 setfl(bflags, bci);
259 gwlarson 3.1 /* run through growing queue */
260     for (bqc = 0; bqc < bqlen; bqc++) {
261 gwlarson 3.2 /* add neighbors until full */
262     for (i = firstneigh(hinp,bqueue[bqc]); i > 0;
263     i = nextneigh()) {
264     if (isset(bflags, i)) /* done already? */
265 gwlarson 3.1 continue;
266 gwlarson 3.2 bqueue[bqlen++] = i; /* add it */
267     bqtotal += bnrays(hinp, i);
268     setfl(bflags, i);
269 gwlarson 3.1 if (bqtotal >= BKBSIZE*1024/sizeof(RAYVAL))
270     break; /* queue full */
271     }
272 gwlarson 3.2 if (i > 0)
273     break;
274 gwlarson 3.1 }
275 gwlarson 3.2 beamdir = hinp->bi; /* sort queue */
276     qsort((char *)bqueue, bqlen, sizeof(*bqueue), bpcmp);
277     /* transfer each beam */
278     for (i = 0; i < bqlen; i++) {
279     bp = hdgetbeam(hinp, bqueue[i]);
280     bcopy((char *)hdbray(bp),
281     (char *)hdnewrays(hout,bqueue[i],bp->nrm),
282     bp->nrm*sizeof(RAYVAL));
283     hdfreebeam(hinp, bqueue[i]);
284     }
285 gwlarson 3.1 hdfreebeam(hout, 0); /* flush output block */
286 gwlarson 3.2 #ifdef DEBUG
287     hdsync(hout, 0);
288     #endif
289 gwlarson 3.1 }
290     /* we're done -- clean up */
291     free((char *)bqueue);
292 gwlarson 3.2 free((char *)bflags);
293 gwlarson 3.1 hddone(hinp);
294     hddone(hout);
295     }
296    
297    
298     eputs(s) /* put error message to stderr */
299     register char *s;
300     {
301     static int midline = 0;
302    
303     if (!*s)
304     return;
305     if (!midline++) { /* prepend line with program name */
306     fputs(progname, stderr);
307     fputs(": ", stderr);
308     }
309     fputs(s, stderr);
310     if (s[strlen(s)-1] == '\n') {
311     fflush(stderr);
312     midline = 0;
313     }
314     }
315    
316    
317     quit(code) /* exit the program gracefully */
318     int code;
319     {
320     hdsync(NULL, 1); /* write out any buffered data */
321     exit(code);
322     }