Back to index

radiance  4R0+20100331
raytrace.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char RCSid[] = "$Id: raytrace.c,v 2.60 2009/12/12 00:03:42 greg Exp $";
00003 #endif
00004 /*
00005  *  raytrace.c - routines for tracing and shading rays.
00006  *
00007  *  External symbols declared in ray.h
00008  */
00009 
00010 #include "copyright.h"
00011 
00012 #include  "ray.h"
00013 #include  "source.h"
00014 #include  "otypes.h"
00015 #include  "otspecial.h"
00016 #include  "random.h"
00017 
00018 #define  MAXCSET     ((MAXSET+1)*2-1)     /* maximum check set size */
00019 
00020 RNUMBER  raynum = 0;        /* next unique ray number */
00021 RNUMBER  nrays = 0;         /* number of calls to localhit */
00022 
00023 static RREAL  Lambfa[5] = {PI, PI, PI, 0.0, 0.0};
00024 OBJREC  Lamb = {
00025        OVOID, MAT_PLASTIC, "Lambertian",
00026        {NULL, Lambfa, 0, 5}, NULL
00027 };                                 /* a Lambertian surface */
00028 
00029 OBJREC  Aftplane;                  /* aft clipping plane object */
00030 
00031 #define  RAYHIT             (-1)          /* return value for intercepted ray */
00032 
00033 static int raymove(FVECT  pos, OBJECT  *cxs, int  dirf, RAY  *r, CUBE  *cu);
00034 static int checkhit(RAY  *r, CUBE  *cu, OBJECT  *cxs);
00035 static void checkset(OBJECT  *os, OBJECT  *cs);
00036 
00037 
00038 extern int
00039 rayorigin(           /* start new ray from old one */
00040        RAY  *r,
00041        int  rt,
00042        const RAY  *ro,
00043        const COLOR rc
00044 )
00045 {
00046        double rw, re;
00047                                           /* assign coefficient/weight */
00048        if (rc == NULL) {
00049               rw = 1.0;
00050               setcolor(r->rcoef, 1., 1., 1.);
00051        } else {
00052               rw = intens(rc);
00053               if (rc != r->rcoef)
00054                      copycolor(r->rcoef, rc);
00055        }
00056        if ((r->parent = ro) == NULL) {           /* primary ray */
00057               r->rlvl = 0;
00058               r->rweight = rw;
00059               r->crtype = r->rtype = rt;
00060               r->rsrc = -1;
00061               r->clipset = NULL;
00062               r->revf = raytrace;
00063               copycolor(r->cext, cextinction);
00064               copycolor(r->albedo, salbedo);
00065               r->gecc = seccg;
00066               r->slights = NULL;
00067        } else {                           /* spawned ray */
00068               if (ro->rot >= FHUGE) {
00069                      memset(r, 0, sizeof(RAY));
00070                      return(-1);          /* illegal continuation */
00071               }
00072               r->rlvl = ro->rlvl;
00073               if (rt & RAYREFL) {
00074                      r->rlvl++;
00075                      r->rsrc = -1;
00076                      r->clipset = ro->clipset;
00077                      r->rmax = 0.0;
00078               } else {
00079                      r->rsrc = ro->rsrc;
00080                      r->clipset = ro->newcset;
00081                      r->rmax = ro->rmax <= FTINY ? 0.0 : ro->rmax - ro->rot;
00082               }
00083               r->revf = ro->revf;
00084               copycolor(r->cext, ro->cext);
00085               copycolor(r->albedo, ro->albedo);
00086               r->gecc = ro->gecc;
00087               r->slights = ro->slights;
00088               r->crtype = ro->crtype | (r->rtype = rt);
00089               VCOPY(r->rorg, ro->rop);
00090               r->rweight = ro->rweight * rw;
00091                                           /* estimate extinction */
00092               re = colval(ro->cext,RED) < colval(ro->cext,GRN) ?
00093                             colval(ro->cext,RED) : colval(ro->cext,GRN);
00094               if (colval(ro->cext,BLU) < re) re = colval(ro->cext,BLU);
00095               re *= ro->rot;
00096               if (re > 0.1) {
00097                      if (re > 92.) {
00098                             r->rweight = 0.0;
00099                      } else {
00100                             r->rweight *= exp(-re);
00101                      }
00102               }
00103        }
00104        rayclear(r);
00105        if (r->rweight <= 0.0)                    /* check for expiration */
00106               return(-1);
00107        if (r->crtype & SHADOW)                   /* shadow commitment */
00108               return(0);
00109        if (maxdepth <= 0 && rc != NULL) { /* Russian roulette */
00110               if (minweight <= 0.0)
00111                      error(USER, "zero ray weight in Russian roulette");
00112               if (maxdepth < 0 && r->rlvl > -maxdepth)
00113                      return(-1);          /* upper reflection limit */
00114               if (r->rweight >= minweight)
00115                      return(0);
00116               if (frandom() > r->rweight/minweight)
00117                      return(-1);
00118               rw = minweight/r->rweight;  /* promote survivor */
00119               scalecolor(r->rcoef, rw);
00120               r->rweight = minweight;
00121               return(0);
00122        }
00123        return(r->rlvl <= abs(maxdepth) && r->rweight >= minweight ? 0 : -1);
00124 }
00125 
00126 
00127 extern void
00128 rayclear(                   /* clear a ray for (re)evaluation */
00129        register RAY  *r
00130 )
00131 {
00132        r->rno = raynum++;
00133        r->newcset = r->clipset;
00134        r->hitf = rayhit;
00135        r->robj = OVOID;
00136        r->ro = NULL;
00137        r->rox = NULL;
00138        r->rt = r->rot = FHUGE;
00139        r->pert[0] = r->pert[1] = r->pert[2] = 0.0;
00140        r->uv[0] = r->uv[1] = 0.0;
00141        setcolor(r->pcol, 1.0, 1.0, 1.0);
00142        setcolor(r->rcol, 0.0, 0.0, 0.0);
00143 }
00144 
00145 
00146 extern void
00147 raytrace(                   /* trace a ray and compute its value */
00148        RAY  *r
00149 )
00150 {
00151        if (localhit(r, &thescene))
00152               raycont(r);          /* hit local surface, evaluate */
00153        else if (r->ro == &Aftplane) {
00154               r->ro = NULL;        /* hit aft clipping plane */
00155               r->rot = FHUGE;
00156        } else if (sourcehit(r))
00157               rayshade(r, r->ro->omod);   /* distant source */
00158 
00159        if (trace != NULL)
00160               (*trace)(r);         /* trace execution */
00161 
00162        rayparticipate(r);          /* for participating medium */
00163 }
00164 
00165 
00166 extern void
00167 raycont(                    /* check for clipped object and continue */
00168        register RAY  *r
00169 )
00170 {
00171        if ((r->clipset != NULL && inset(r->clipset, r->ro->omod)) ||
00172                      !rayshade(r, r->ro->omod))
00173               raytrans(r);
00174 }
00175 
00176 
00177 extern void
00178 raytrans(                   /* transmit ray as is */
00179        register RAY  *r
00180 )
00181 {
00182        RAY  tr;
00183 
00184        if (rayorigin(&tr, TRANS, r, NULL) == 0) {
00185               VCOPY(tr.rdir, r->rdir);
00186               rayvalue(&tr);
00187               copycolor(r->rcol, tr.rcol);
00188               r->rt = r->rot + tr.rt;
00189        }
00190 }
00191 
00192 
00193 extern int
00194 rayshade(            /* shade ray r with material mod */
00195        register RAY  *r,
00196        int  mod
00197 )
00198 {
00199        register OBJREC  *m;
00200 
00201        r->rt = r->rot;                    /* set effective ray length */
00202        for ( ; mod != OVOID; mod = m->omod) {
00203               m = objptr(mod);
00204               /****** unnecessary test since modifier() is always called
00205               if (!ismodifier(m->otype)) {
00206                      sprintf(errmsg, "illegal modifier \"%s\"", m->oname);
00207                      error(USER, errmsg);
00208               }
00209               ******/
00210                                    /* hack for irradiance calculation */
00211               if (do_irrad && !(r->crtype & ~(PRIMARY|TRANS)) &&
00212                             m->otype != MAT_CLIP &&
00213                             (ofun[m->otype].flags & (T_M|T_X))) {
00214                      if (irr_ignore(m->otype)) {
00215                             raytrans(r);
00216                             return(1);
00217                      }
00218                      if (!islight(m->otype))
00219                             m = &Lamb;
00220               }
00221               if ((*ofun[m->otype].funp)(m, r))
00222                      return(1);    /* materials call raytexture() */
00223        }
00224        return(0);                  /* no material! */
00225 }
00226 
00227 
00228 extern void
00229 rayparticipate(                    /* compute ray medium participation */
00230        register RAY  *r
00231 )
00232 {
00233        COLOR  ce, ca;
00234        double re, ge, be;
00235 
00236        if (intens(r->cext) <= 1./FHUGE)
00237               return;                            /* no medium */
00238        re = r->rot*colval(r->cext,RED);
00239        ge = r->rot*colval(r->cext,GRN);
00240        be = r->rot*colval(r->cext,BLU);
00241        if (r->crtype & SHADOW) {          /* no scattering for sources */
00242               re *= 1. - colval(r->albedo,RED);
00243               ge *= 1. - colval(r->albedo,GRN);
00244               be *= 1. - colval(r->albedo,BLU);
00245        }
00246        setcolor(ce,  re<=FTINY ? 1. : re>92. ? 0. : exp(-re),
00247                      ge<=FTINY ? 1. : ge>92. ? 0. : exp(-ge),
00248                      be<=FTINY ? 1. : be>92. ? 0. : exp(-be));
00249        multcolor(r->rcol, ce);                   /* path extinction */
00250        if (r->crtype & SHADOW || intens(r->albedo) <= FTINY)
00251               return;                            /* no scattering */
00252        setcolor(ca,
00253               colval(r->albedo,RED)*colval(ambval,RED)*(1.-colval(ce,RED)),
00254               colval(r->albedo,GRN)*colval(ambval,GRN)*(1.-colval(ce,GRN)),
00255               colval(r->albedo,BLU)*colval(ambval,BLU)*(1.-colval(ce,BLU)));
00256        addcolor(r->rcol, ca);                    /* ambient in scattering */
00257        srcscatter(r);                            /* source in scattering */
00258 }
00259 
00260 
00261 extern void
00262 raytexture(                 /* get material modifiers */
00263        RAY  *r,
00264        OBJECT  mod
00265 )
00266 {
00267        register OBJREC  *m;
00268                                    /* execute textures and patterns */
00269        for ( ; mod != OVOID; mod = m->omod) {
00270               m = objptr(mod);
00271               /****** unnecessary test since modifier() is always called
00272               if (!ismodifier(m->otype)) {
00273                      sprintf(errmsg, "illegal modifier \"%s\"", m->oname);
00274                      error(USER, errmsg);
00275               }
00276               ******/
00277               if ((*ofun[m->otype].funp)(m, r)) {
00278                      sprintf(errmsg, "conflicting material \"%s\"",
00279                                    m->oname);
00280                      objerror(r->ro, USER, errmsg);
00281               }
00282        }
00283 }
00284 
00285 
00286 extern int
00287 raymixture(          /* mix modifiers */
00288        register RAY  *r,
00289        OBJECT  fore,
00290        OBJECT  back,
00291        double  coef
00292 )
00293 {
00294        RAY  fr, br;
00295        int  foremat, backmat;
00296        register int  i;
00297                                    /* bound coefficient */
00298        if (coef > 1.0)
00299               coef = 1.0;
00300        else if (coef < 0.0)
00301               coef = 0.0;
00302                                    /* compute foreground and background */
00303        foremat = backmat = 0;
00304                                    /* foreground */
00305        fr = *r;
00306        if (coef > FTINY) {
00307               fr.rweight *= coef;
00308               scalecolor(fr.rcoef, coef);
00309               foremat = rayshade(&fr, fore);
00310        }
00311                                    /* background */
00312        br = *r;
00313        if (coef < 1.0-FTINY) {
00314               br.rweight *= 1.0-coef;
00315               scalecolor(br.rcoef, 1.0-coef);
00316               backmat = rayshade(&br, back);
00317        }
00318                                    /* check for transparency */
00319        if (backmat ^ foremat) {
00320               if (backmat && coef > FTINY)
00321                      raytrans(&fr);
00322               else if (foremat && coef < 1.0-FTINY)
00323                      raytrans(&br);
00324        }
00325                                    /* mix perturbations */
00326        for (i = 0; i < 3; i++)
00327               r->pert[i] = coef*fr.pert[i] + (1.0-coef)*br.pert[i];
00328                                    /* mix pattern colors */
00329        scalecolor(fr.pcol, coef);
00330        scalecolor(br.pcol, 1.0-coef);
00331        copycolor(r->pcol, fr.pcol);
00332        addcolor(r->pcol, br.pcol);
00333                                    /* return value tells if material */
00334        if (!foremat & !backmat)
00335               return(0);
00336                                    /* mix returned ray values */
00337        scalecolor(fr.rcol, coef);
00338        scalecolor(br.rcol, 1.0-coef);
00339        copycolor(r->rcol, fr.rcol);
00340        addcolor(r->rcol, br.rcol);
00341        r->rt = bright(fr.rcol) > bright(br.rcol) ? fr.rt : br.rt;
00342        return(1);
00343 }
00344 
00345 
00346 extern double
00347 raydist(             /* compute (cumulative) ray distance */
00348        register const RAY  *r,
00349        register int  flags
00350 )
00351 {
00352        double  sum = 0.0;
00353 
00354        while (r != NULL && r->crtype&flags) {
00355               sum += r->rot;
00356               r = r->parent;
00357        }
00358        return(sum);
00359 }
00360 
00361 
00362 extern void
00363 raycontrib(          /* compute (cumulative) ray contribution */
00364        double  rc[3],
00365        const RAY  *r,
00366        int  flags
00367 )
00368 {
00369        double eext[3];
00370        int    i;
00371 
00372        eext[0] = eext[1] = eext[2] = 0.;
00373        rc[0] = rc[1] = rc[2] = 1.;
00374 
00375        while (r != NULL && r->crtype&flags) {
00376               for (i = 3; i--; ) {
00377                      rc[i] *= colval(r->rcoef,i);
00378                      eext[i] += r->rot * colval(r->cext,i);
00379               }
00380               r = r->parent;
00381        }
00382        for (i = 3; i--; )
00383               rc[i] *= (eext[i] <= FTINY) ? 1. :
00384                             (eext[i] > 92.) ? 0. : exp(-eext[i]);
00385 }
00386 
00387 
00388 extern double
00389 raynormal(           /* compute perturbed normal for ray */
00390        FVECT  norm,
00391        register RAY  *r
00392 )
00393 {
00394        double  newdot;
00395        register int  i;
00396 
00397        /*     The perturbation is added to the surface normal to obtain
00398         *  the new normal.  If the new normal would affect the surface
00399         *  orientation wrt. the ray, a correction is made.  The method is
00400         *  still fraught with problems since reflected rays and similar
00401         *  directions calculated from the surface normal may spawn rays behind
00402         *  the surface.  The only solution is to curb textures at high
00403         *  incidence (namely, keep DOT(rdir,pert) < Rdot).
00404         */
00405 
00406        for (i = 0; i < 3; i++)
00407               norm[i] = r->ron[i] + r->pert[i];
00408 
00409        if (normalize(norm) == 0.0) {
00410               objerror(r->ro, WARNING, "illegal normal perturbation");
00411               VCOPY(norm, r->ron);
00412               return(r->rod);
00413        }
00414        newdot = -DOT(norm, r->rdir);
00415        if ((newdot > 0.0) != (r->rod > 0.0)) {          /* fix orientation */
00416               for (i = 0; i < 3; i++)
00417                      norm[i] += 2.0*newdot*r->rdir[i];
00418               newdot = -newdot;
00419        }
00420        return(newdot);
00421 }
00422 
00423 
00424 extern void
00425 newrayxf(                   /* get new tranformation matrix for ray */
00426        RAY  *r
00427 )
00428 {
00429        static struct xfn {
00430               struct xfn  *next;
00431               FULLXF  xf;
00432        }  xfseed = { &xfseed }, *xflast = &xfseed;
00433        register struct xfn  *xp;
00434        register const RAY  *rp;
00435 
00436        /*
00437         * Search for transform in circular list that
00438         * has no associated ray in the tree.
00439         */
00440        xp = xflast;
00441        for (rp = r->parent; rp != NULL; rp = rp->parent)
00442               if (rp->rox == &xp->xf) {          /* xp in use */
00443                      xp = xp->next;                     /* move to next */
00444                      if (xp == xflast) {         /* need new one */
00445                             xp = (struct xfn *)malloc(sizeof(struct xfn));
00446                             if (xp == NULL)
00447                                    error(SYSTEM,
00448                                           "out of memory in newrayxf");
00449                                                  /* insert in list */
00450                             xp->next = xflast->next;
00451                             xflast->next = xp;
00452                             break;               /* we're done */
00453                      }
00454                      rp = r;                     /* start check over */
00455               }
00456                                    /* got it */
00457        r->rox = &xp->xf;
00458        xflast = xp;
00459 }
00460 
00461 
00462 extern void
00463 flipsurface(                /* reverse surface orientation */
00464        register RAY  *r
00465 )
00466 {
00467        r->rod = -r->rod;
00468        r->ron[0] = -r->ron[0];
00469        r->ron[1] = -r->ron[1];
00470        r->ron[2] = -r->ron[2];
00471        r->pert[0] = -r->pert[0];
00472        r->pert[1] = -r->pert[1];
00473        r->pert[2] = -r->pert[2];
00474 }
00475 
00476 
00477 extern void
00478 rayhit(                     /* standard ray hit test */
00479        OBJECT  *oset,
00480        RAY  *r
00481 )
00482 {
00483        OBJREC  *o;
00484        int    i;
00485 
00486        for (i = oset[0]; i > 0; i--) {
00487               o = objptr(oset[i]);
00488               if ((*ofun[o->otype].funp)(o, r))
00489                      r->robj = oset[i];
00490        }
00491 }
00492 
00493 
00494 extern int
00495 localhit(            /* check for hit in the octree */
00496        register RAY  *r,
00497        register CUBE  *scene
00498 )
00499 {
00500        OBJECT  cxset[MAXCSET+1];   /* set of checked objects */
00501        FVECT  curpos;                     /* current cube position */
00502        int  sflags;                /* sign flags */
00503        double  t, dt;
00504        register int  i;
00505 
00506        nrays++;                    /* increment trace counter */
00507        sflags = 0;
00508        for (i = 0; i < 3; i++) {
00509               curpos[i] = r->rorg[i];
00510               if (r->rdir[i] > 1e-7)
00511                      sflags |= 1 << i;
00512               else if (r->rdir[i] < -1e-7)
00513                      sflags |= 0x10 << i;
00514        }
00515        if (sflags == 0)
00516               error(CONSISTENCY, "zero ray direction in localhit");
00517                                    /* start off assuming nothing hit */
00518        if (r->rmax > FTINY) {             /* except aft plane if one */
00519               r->ro = &Aftplane;
00520               r->rot = r->rmax;
00521               for (i = 0; i < 3; i++)
00522                      r->rop[i] = r->rorg[i] + r->rot*r->rdir[i];
00523        }
00524                                    /* find global cube entrance point */
00525        t = 0.0;
00526        if (!incube(scene, curpos)) {
00527                                    /* find distance to entry */
00528               for (i = 0; i < 3; i++) {
00529                                    /* plane in our direction */
00530                      if (sflags & 1<<i)
00531                             dt = scene->cuorg[i];
00532                      else if (sflags & 0x10<<i)
00533                             dt = scene->cuorg[i] + scene->cusize;
00534                      else
00535                             continue;
00536                                    /* distance to the plane */
00537                      dt = (dt - r->rorg[i])/r->rdir[i];
00538                      if (dt > t)
00539                             t = dt;       /* farthest face is the one */
00540               }
00541               t += FTINY;          /* fudge to get inside cube */
00542               if (t >= r->rot)     /* clipped already */
00543                      return(0);
00544                                    /* advance position */
00545               for (i = 0; i < 3; i++)
00546                      curpos[i] += r->rdir[i]*t;
00547 
00548               if (!incube(scene, curpos)) /* non-intersecting ray */
00549                      return(0);
00550        }
00551        cxset[0] = 0;
00552        raymove(curpos, cxset, sflags, r, scene);
00553        return((r->ro != NULL) & (r->ro != &Aftplane));
00554 }
00555 
00556 
00557 static int
00558 raymove(             /* check for hit as we move */
00559        FVECT  pos,                 /* current position, modified herein */
00560        OBJECT  *cxs,               /* checked objects, modified by checkhit */
00561        int  dirf,                  /* direction indicators to speed tests */
00562        register RAY  *r,
00563        register CUBE  *cu
00564 )
00565 {
00566        int  ax;
00567        double  dt, t;
00568 
00569        if (istree(cu->cutree)) {          /* recurse on subcubes */
00570               CUBE  cukid;
00571               register int  br, sgn;
00572 
00573               cukid.cusize = cu->cusize * 0.5;   /* find subcube */
00574               VCOPY(cukid.cuorg, cu->cuorg);
00575               br = 0;
00576               if (pos[0] >= cukid.cuorg[0]+cukid.cusize) {
00577                      cukid.cuorg[0] += cukid.cusize;
00578                      br |= 1;
00579               }
00580               if (pos[1] >= cukid.cuorg[1]+cukid.cusize) {
00581                      cukid.cuorg[1] += cukid.cusize;
00582                      br |= 2;
00583               }
00584               if (pos[2] >= cukid.cuorg[2]+cukid.cusize) {
00585                      cukid.cuorg[2] += cukid.cusize;
00586                      br |= 4;
00587               }
00588               for ( ; ; ) {
00589                      cukid.cutree = octkid(cu->cutree, br);
00590                      if ((ax = raymove(pos,cxs,dirf,r,&cukid)) == RAYHIT)
00591                             return(RAYHIT);
00592                      sgn = 1 << ax;
00593                      if (sgn & dirf)                    /* positive axis? */
00594                             if (sgn & br)
00595                                    return(ax);   /* overflow */
00596                             else {
00597                                    cukid.cuorg[ax] += cukid.cusize;
00598                                    br |= sgn;
00599                             }
00600                      else
00601                             if (sgn & br) {
00602                                    cukid.cuorg[ax] -= cukid.cusize;
00603                                    br &= ~sgn;
00604                             } else
00605                                    return(ax);   /* underflow */
00606               }
00607               /*NOTREACHED*/
00608        }
00609        if (isfull(cu->cutree)) {
00610               if (checkhit(r, cu, cxs))
00611                      return(RAYHIT);
00612        } else if (r->ro == &Aftplane && incube(cu, r->rop))
00613               return(RAYHIT);
00614                                    /* advance to next cube */
00615        if (dirf&0x11) {
00616               dt = dirf&1 ? cu->cuorg[0] + cu->cusize : cu->cuorg[0];
00617               t = (dt - pos[0])/r->rdir[0];
00618               ax = 0;
00619        } else
00620               t = FHUGE;
00621        if (dirf&0x22) {
00622               dt = dirf&2 ? cu->cuorg[1] + cu->cusize : cu->cuorg[1];
00623               dt = (dt - pos[1])/r->rdir[1];
00624               if (dt < t) {
00625                      t = dt;
00626                      ax = 1;
00627               }
00628        }
00629        if (dirf&0x44) {
00630               dt = dirf&4 ? cu->cuorg[2] + cu->cusize : cu->cuorg[2];
00631               dt = (dt - pos[2])/r->rdir[2];
00632               if (dt < t) {
00633                      t = dt;
00634                      ax = 2;
00635               }
00636        }
00637        pos[0] += r->rdir[0]*t;
00638        pos[1] += r->rdir[1]*t;
00639        pos[2] += r->rdir[2]*t;
00640        return(ax);
00641 }
00642 
00643 
00644 static int
00645 checkhit(            /* check for hit in full cube */
00646        register RAY  *r,
00647        CUBE  *cu,
00648        OBJECT  *cxs
00649 )
00650 {
00651        OBJECT  oset[MAXSET+1];
00652 
00653        objset(oset, cu->cutree);
00654        checkset(oset, cxs);               /* avoid double-checking */
00655 
00656        (*r->hitf)(oset, r);               /* test for hit in set */
00657 
00658        if (r->robj == OVOID)
00659               return(0);                  /* no scores yet */
00660 
00661        return(incube(cu, r->rop));        /* hit OK if in current cube */
00662 }
00663 
00664 
00665 static void
00666 checkset(            /* modify checked set and set to check */
00667        register OBJECT  *os,                     /* os' = os - cs */
00668        register OBJECT  *cs               /* cs' = cs + os */
00669 )
00670 {
00671        OBJECT  cset[MAXCSET+MAXSET+1];
00672        register int  i, j;
00673        int  k;
00674                                    /* copy os in place, cset <- cs */
00675        cset[0] = 0;
00676        k = 0;
00677        for (i = j = 1; i <= os[0]; i++) {
00678               while (j <= cs[0] && cs[j] < os[i])
00679                      cset[++cset[0]] = cs[j++];
00680               if (j > cs[0] || os[i] != cs[j]) { /* object to check */
00681                      os[++k] = os[i];
00682                      cset[++cset[0]] = os[i];
00683               }
00684        }
00685        if (!(os[0] = k))           /* new "to check" set size */
00686               return;                     /* special case */
00687        while (j <= cs[0])          /* get the rest of cs */
00688               cset[++cset[0]] = cs[j++];
00689        if (cset[0] > MAXCSET)             /* truncate "checked" set if nec. */
00690               cset[0] = MAXCSET;
00691        /* setcopy(cs, cset); */    /* copy cset back to cs */
00692        os = cset;
00693        for (i = os[0]; i-- >= 0; )
00694               *cs++ = *os++;
00695 }