Back to index

radiance  4R0+20100331
clumpbeams.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: clumpbeams.c,v 3.8 2004/01/01 11:21:55 schorsch Exp $";
00003 #endif
00004 /*
00005  * Bundle holodeck beams together into clumps.
00006  */
00007 
00008 #include "holo.h"
00009 
00010 #define flgop(p,i,op)              ((p)[(i)>>5] op (1L<<((i)&0x1f)))
00011 #define issetfl(p,i)        flgop(p,i,&)
00012 #define setfl(p,i)          flgop(p,i,|=)
00013 #define clrfl(p,i)          flgop(p,i,&=~)
00014 
00015 static int    bneighlist[9*9-1];
00016 static int    bneighrem;
00017 
00018 #define nextneigh()  (bneighrem<=0 ? 0 : bneighlist[--bneighrem])
00019 
00020 static void gcshifti(GCOORD *gc, int ia, int di, HOLO *hp);
00021 static void mkneighgrid(GCOORD     ng[3*3], HOLO *hp, GCOORD   *gc);
00022 static int firstneigh(HOLO  *hp, int      b);
00023 
00024 
00025 
00026 static void
00027 gcshifti(     /* shift cell row or column */
00028        register GCOORD      *gc,
00029        int    ia,
00030        int    di,
00031        register HOLO *hp
00032 )
00033 {
00034        int    nw;
00035 
00036        if (di > 0) {
00037               if (++gc->i[ia] >= hp->grid[((gc->w>>1)+1+ia)%3]) {
00038                      nw = ((gc->w&~1) + (ia<<1) + 3) % 6;
00039                      gc->i[ia] = gc->i[1-ia];
00040                      gc->i[1-ia] = gc->w&1 ? hp->grid[((nw>>1)+2-ia)%3]-1 : 0;
00041                      gc->w = nw;
00042               }
00043        } else if (di < 0) {
00044               if (--gc->i[ia] < 0) {
00045                      nw = ((gc->w&~1) + (ia<<1) + 2) % 6;
00046                      gc->i[ia] = gc->i[1-ia];
00047                      gc->i[1-ia] = gc->w&1 ? hp->grid[((nw>>1)+2-ia)%3]-1 : 0;
00048                      gc->w = nw;
00049               }
00050        }
00051 }
00052 
00053 
00054 static void
00055 mkneighgrid(         /* compute neighborhood for grid cell */
00056        GCOORD ng[3*3],
00057        HOLO   *hp,
00058        GCOORD *gc
00059 )
00060 {
00061        GCOORD gci0;
00062        register int  i, j;
00063 
00064        for (i = 3; i--; ) {
00065               gci0 = *gc;
00066               gcshifti(&gci0, 0, i-1, hp);
00067               for (j = 3; j--; ) {
00068                      *(ng+(3*i+j)) = gci0;
00069                      gcshifti(ng+(3*i+j), gci0.w==gc->w, j-1, hp);
00070               }
00071        }
00072 }
00073 
00074 
00075 static int
00076 firstneigh(          /* initialize neighbor list and return first */
00077        HOLO   *hp,
00078        int    b
00079 )
00080 {
00081        GCOORD wg0[9], wg1[9], bgc[2];
00082        int    i, j;
00083 
00084        hdbcoord(bgc, hp, b);
00085        mkneighgrid(wg0, hp, bgc);
00086        mkneighgrid(wg1, hp, bgc+1);
00087        bneighrem = 0;
00088        for (i = 9; i--; )
00089               for (j = 9; j--; ) {
00090                      if ((i == 4) & (j == 4))    /* don't copy starting beam */
00091                             continue;
00092                      if (wg0[i].w == wg1[j].w)
00093                             continue;
00094                      *bgc = *(wg0+i);
00095                      *(bgc+1) = *(wg1+j);
00096                      bneighlist[bneighrem++] = hdbindex(hp, bgc);
00097 #ifdef DEBUG
00098                      if (bneighlist[bneighrem-1] <= 0)
00099                             error(CONSISTENCY, "bad beam in firstneigh");
00100 #endif
00101               }
00102        return(nextneigh());
00103 }
00104 
00105 
00106 extern void
00107 clumpbeams(   /* clump beams from hinp */
00108        register HOLO *hp,
00109        int    maxcnt,
00110        int    maxsiz,
00111        int    (*cf)(HOLO *hp, int *bqueue, int bqlen)
00112 )
00113 {
00114        static short  primes[] = {9431,6803,4177,2659,1609,887,587,251,47,1};
00115        uint32 *bflags;
00116        int    *bqueue;
00117        int    bqlen;
00118        int32  bqtotal;
00119        int    bc, bci, bqc, myprime;
00120        register int  i;
00121                                    /* get clump size */
00122        if (maxcnt <= 1)
00123               maxcnt = nbeams(hp);
00124        maxsiz /= sizeof(RAYVAL);
00125                                    /* allocate beam queue */
00126        bqueue = (int *)malloc(maxcnt*sizeof(int));
00127        bflags = (uint32 *)calloc((nbeams(hp)>>5)+1,
00128                      sizeof(uint32));
00129        if ((bqueue == NULL) | (bflags == NULL))
00130               error(SYSTEM, "out of memory in clumpbeams");
00131                                    /* mark empty beams as done */
00132        for (i = nbeams(hp); i > 0; i--)
00133               if (!bnrays(hp, i))
00134                      setfl(bflags, i);
00135                                    /* pick a good prime step size */
00136        for (i = 0; primes[i]<<5 >= nbeams(hp); i++)
00137               ;
00138        while ((myprime = primes[i++]) > 1)
00139               if (nbeams(hp) % myprime)
00140                      break;
00141                                    /* add each input beam and neighbors */
00142        for (bc = bci = nbeams(hp); bc > 0; bc--,
00143                      bci += bci>myprime ? -myprime : nbeams(hp)-myprime) {
00144               if (issetfl(bflags, bci))
00145                      continue;
00146               bqueue[0] = bci;            /* initialize queue */
00147               bqlen = 1;
00148               bqtotal = bnrays(hp, i);
00149               setfl(bflags, bci);
00150                                           /* run through growing queue */
00151               for (bqc = 0; bqc < bqlen; bqc++) {
00152                                           /* add neighbors until full */
00153                      for (i = firstneigh(hp,bqueue[bqc]); i > 0;
00154                                    i = nextneigh()) {
00155                             if (issetfl(bflags, i))     /* done already? */
00156                                    continue;
00157                             bqueue[bqlen++] = i; /* add it */
00158                             bqtotal += bnrays(hp, i);
00159                             setfl(bflags, i);
00160                             if (bqlen >= maxcnt ||
00161                                           (maxsiz && bqtotal >= maxsiz))
00162                                    break;        /* queue full */
00163                      }
00164                      if (i > 0)
00165                             break;
00166               }
00167               (*cf)(hp, bqueue, bqlen);   /* transfer clump */
00168        }
00169                                    /* all done; clean up */
00170        free((void *)bqueue);
00171        free((void *)bflags);
00172 }