zbufsort.c

很经典的电磁计算（电容计算）软件 MIT90年代开发

  /* a larger alpha => greater distance from view plane */
  /* => fac is deeper than facref if alpha_fac > alpha_facref */
  /*     (ie if alpha_facref - alpha_fac < 0) */
  /* => facref is deeper than fac if alpha_facref > alpha_fac */
  /* => inconclusive if equal */
  if(diff_is_negative(alpha_facref, alpha_fac, 1.0)) return(TRUE);
  else if(diff_is_negative(alpha_fac, alpha_facref, 1.0)) return(REVERSE);
  else {
    fprintf(stderr, 
	    "\nis1stFaceDeeper: Warning, face ordering test failure\n");
    fprintf(stderr, "  alpha_fac, face %d = %g alpha_facref, face %d = %g\n", 
	    fac->index, alpha_fac, facref->index, alpha_facref);
    dump_face(stderr, fac);
    fprintf(stderr, " Projected corners\n");
    dumpCorners(stderr, cproj[0], fac->numsides, 3);
    dump_face(stderr, facref);
    fprintf(stderr, " Projected corners\n");
    dumpCorners(stderr, cproj[1], facref->numsides, 3);

    return(FALSE);		/* inconclusive test */
  }

}
#endif

/*
  returns TRUE if bounding box of facref and fac (proj to facref's plane) insct
*/
int isThereBoxOverlap(fac, facref, view)
face *fac, *facref;
double *view;
{
  int i, j, olap[2];
  double cproj[MAXSIDES][3];	/* corners of fac in facref's plane */
  double cref[MAXSIDES][3];	/* corners of facref in facref plane */
  double alpha[MAXSIDES];	/* 1 => view point 0 => corner */
  double minref[2], maxref[2];	/* bounding box coordinates */
  double minfac[2], maxfac[2];
  double x[3], y[3];		/* coordinates of x and y in facref plane */
  double dot(), temp, tvec[3], tvec1[3], margin, ovrlapmgn = 0.0;

  /* figure projections of fac's corners back to facref's plane rel to view */
  for(i = 0; i < fac->numsides; i++) {
    for(j = 0; j < 3; j++) tvec[j] = view[j] - fac->c[i][j]; /* get v-c */
    temp = dot(facref->normal, tvec); /* get n.(v-c) */
    margin = sqrt(dot(tvec, tvec))*MARGIN;
    ovrlapmgn = MAX(margin, ovrlapmgn);	/* used below */
    /* test fails if v-c is perpendicular to n */
    if(temp > -margin && temp < margin) return(FALSE);
    /* get alpha as in c + alpha*(v-c) = c' */
    alpha[i] = (facref->rhs - dot(facref->normal, fac->c[i]))/temp;
    /* if(alpha[i] < -margin || alpha[i] > 1.0+margin) {
      fprintf(stderr, "\nisThereBoxOverlap: big X failure, alpha = %g\n",
	      alpha[i]);
      exit(1);
    } */
    for(j = 0; j < 3; j++)	/* get c' */
	cproj[i][j] = (1.0-alpha[i])*fac->c[i][j]+alpha[i]*view[j];
  }

  /* figure x and y coordinates in facref plane (normal is z coordinate) */
  /* x = c0-c1 always */
  for(j = 0; j < 3; j++) x[j] = facref->c[0][j] - facref->c[1][j];
  temp = sqrt(dot(x, x));
  for(j = 0; j < 3; j++) x[j] /= temp; /* normalize */
  /* y = zXx */
  crossProd(y, facref->normal, x);
  
  /* project all fac corner projections onto new x and y coordinates
     - facref->c[0] plays the role of origin */
  for(i = 0; i < fac->numsides; i++) {
    for(j = 0; j < 3; j++) tvec1[j] = cproj[i][j] - facref->c[0][j];
    tvec[0] = dot(x, tvec1);	/* get weight in x direction */
    tvec[1] = dot(y, tvec1);	/* get weight in y direction */
    tvec[2] = 0.0;		/* all z weights must = facref->rhs */
    for(j = 0; j < 3; j++) cproj[i][j] = tvec[j]; /* xfer */
  }
  for(j = 0; j < 3; j++) cref[0][j] = 0.0;
  for(i = 1; i < facref->numsides; i++) {
    for(j = 0; j < 3; j++) tvec1[j] = facref->c[i][j] - facref->c[0][j];
    cref[i][0] = dot(x, tvec1);	/* get weight in x direction */
    cref[i][1] = dot(y, tvec1);	/* get weight in y direction */
    cref[i][2] = 0.0;		/* all z weights must = facref->rhs */
  }

  /* figure bounding boxes in new coordinates */
  minfac[0] = maxfac[0] = cproj[0][0];
  minfac[1] = maxfac[1] = cproj[0][1];
  for(i = 1; i < fac->numsides; i++) { /* find max, min for fac */
    for(j = 0; j < 2; j++) {
      minfac[j] = MIN(minfac[j], cproj[i][j]);
      maxfac[j] = MAX(maxfac[j], cproj[i][j]);
    }
  }
  minref[0] = maxref[0] = cref[0][0];
  minref[1] = maxref[1] = cref[0][1];
  for(i = 1; i < facref->numsides; i++) { /* find max, min for facref */
    for(j = 0; j < 2; j++) {
      minref[j] = MIN(minref[j], cref[i][j]);
      maxref[j] = MAX(maxref[j], cref[i][j]);
    }
  }

  /* check for overlap - call things overlaped even if not to be sure
     - no overlap if either no overlap in x or none in y */
  /* check for x overlap */
  olap[0] = olap[1] = FALSE;
  for(j = 0; j < 2; j++) {
    if((minref[j]-ovrlapmgn < minfac[j] && minfac[j] < maxref[j]+ovrlapmgn) ||
       (minref[j]-ovrlapmgn < maxfac[j] && maxfac[j] < maxref[j]+ovrlapmgn) ||
       (minfac[j]-ovrlapmgn < minref[j] && minref[j] < maxfac[j]+ovrlapmgn) ||
       (minfac[j]-ovrlapmgn < maxref[j] && maxref[j] < maxfac[j]+ovrlapmgn))
	olap[j] = TRUE;
  }
  if(olap[0] == FALSE || olap[1] == FALSE) return(FALSE);
  else return(TRUE);
}

#if 1 == 0
/*
  returns TRUE if 1st face is deeper than second - uses Don's 3d X method
  - returns FALSE in don't care situations (1st need not be ordered before 2nd)
*/
int is1stFaceDeeper(first, second, view)
face *first, *second;
double *view;
{
  int fst, snd, ret = SAME;

  /* figure on which side of each face plane lies the other face */
  fst = whichSide(first, second); /* first relative to plane of second */
  snd = whichSide(second, first); /* second relative to plane of first */

  /* figure returned value */
  if((fst == NEG && snd == NEG) || (fst == POS && snd == POS) ||
      (fst == SAME && snd == SAME))
      ret = FALSE;		/* order arbitrary */
  else if((fst == NEG && snd == POS)) ret = TRUE;
  else if((fst == POS && snd == NEG)) ret = FALSE;
  else if(fst == NEG && snd == SPLIT) ret = TRUE;
  else if(fst == POS && snd == SPLIT) ret = FALSE;
  else if(fst == SPLIT && snd == POS) ret = TRUE;
  else if(fst == SPLIT && snd == NEG) ret = FALSE;
  /* next four added 22 Aug 91 w/o really thinking
     - catches false SAME's ? */
/*  else if(fst == NEG && snd == SAME) ret = TRUE;
  else if(fst == POS && snd == SAME) ret = FALSE;
  else if(fst == SAME && snd == POS) ret = TRUE;
  else if(fst == SAME && snd == NEG) ret = FALSE; */
  else if(snd == SPLIT && fst == SPLIT) {
    return(isThereProjOverlap(first, second, view));
#if 1 == 0
    /* if both are split, see if they interact at all */
    if(isThereProjOverlap(first, second, view) == FALSE) ret = FALSE;
    else {
      fprintf(stderr, 
	      "\nis1stFaceDeeper: unresolvable both faces split case\n");
      ret = TRUE;		/* don't order--try to order w/other check
				   sometimes works, sometimes not */
    }
#endif
  }
  if(ret == SAME) {
    fprintf(stderr,"\nis1stFaceDeeper: bad fst = %d and snd = %d\n", fst, snd);
    exit(1);
  }
  else if(ret == TRUE) {	/* make sure 2nd could overlap 1st */
    if(isThereProjOverlap(first, second, view) == FALSE) ret = FALSE;
  }
  return(ret);
}
/*
#else

int is1stFaceDeeper(first, second, view)
face *first, *second;
double *view;
{
  return(isThereProjOverlap(first, second, view));
}
*/
#endif

/*
  recursive guts of below
*/
int chkCycle(fac, ref, fp)
face *fac, *ref;
FILE *fp;
{
  int b, i;

  if(fac->mark == TRUE) return(FALSE);

  fac->mark = TRUE;

  if(fac->numbehind == 0) return(FALSE);
  else {
    for(b = 0; b < fac->numbehind; b++) {
      /*fprintf(fp, " %d (%d)", (fac->behind)[b]->depth,
	      (fac->behind)[b]->index);
      if(b % 5 == 0 && b != 0) fprintf(fp, "\n");*/
      if(fac->behind[b] == ref) return(TRUE);
      else if(chkCycle(fac->behind[b], ref, fp) == TRUE) return(TRUE);
    }
/*    if((i-1) % 5 != 0 || i == 1) fprintf(fp, "\n"); */
  }
  return(FALSE);
}

/*
  checks for cycles in the depth graph - BROKEN(?)
*/
void dumpCycles(faces, numfaces, file)
face **faces;
int numfaces;
FILE *file;
{
  int i, f, j, b, *cycled, *cyclei, numcycle, cycle = FALSE;
  face *fp;

  /* for each face, chase behind pointers until a leaf or same face is found */
  /*fprintf(file, "\nRecursive behind lists\n");*/
  for(f = 0; f < numfaces; f++) {
/*    fprintf(file, "%d (%d):", faces[f]->depth, faces[f]->index);*/
    for(j = 0; j < numfaces; j++) faces[j]->mark = FALSE;
    for(b = 0; b < faces[f]->numbehind; b++) {
      if(chkCycle(faces[f]->behind[b], faces[f], file) == TRUE) {
	cycle = TRUE;
	break;
      }
    }
    if(cycle == TRUE) break;
  }
  if(cycle == FALSE) fprintf(file, "Adjacency graph has no cycles\n");
  else fprintf(file, "Adjacency graph has cycles\n");
  for(j = 0; j < numfaces; j++) faces[j]->mark = FALSE;
}

    
/*
  recursively sets depths of faces
*/
void setDepth(fac)
face *fac;
{
  int i;

  /* mark so this face won't be renumbered */
  fac->mark = TRUE;

  if(fac->index == 193) {
    i = i + 0;
  }

  /* do adjacents if needed */
  for(i = 0; i < fac->numbehind; i++) {
    if((fac->behind[i])->mark == FALSE) setDepth(fac->behind[i]);
  }

  if(fac->index == 131 || fac->index == 193) {
    i = i + 0;
  }

  /* set depth, update counter */
  fac->depth = cnt--;
}

/*
  does a topological sorting of the faces using graph setup by getAdjGraph()
  - returns a new set of pointers with deepest (1st to print) face first
*/
face **depthSortFaces(faces, numfaces)
face **faces;
int numfaces;
{
  int f, i, facefound;
  face **rfaces;

  /* make sure all marks are cleared */
  for(f = 0; f < numfaces; f++) faces[f]->mark = FALSE;

  /* set depths recursively - zero is deepest (first to be rendered) */
  for(f = 0, cnt = numfaces-1; f < numfaces; f++) {
    if(faces[f]->mark == FALSE) setDepth(faces[f]);
  }

  /* make the new set of pointers */
  CALLOC(rfaces, numfaces, face *, ON, AMSC);
  for(f = 0; f < numfaces; f++) {
    for(i = 0, facefound = FALSE; i < numfaces; i++) {
      if(faces[i]->depth == f) {
	rfaces[f] = faces[i];
	facefound = TRUE;
	break;
      }
    }
    if(facefound == FALSE) {
      fprintf(stderr, "\ndepthSortFaces: can't find depth %d face\n", f);
      exit(1);
    }
  }
  return(rfaces);
}

/*
  sets up adjacency graph pointers in faces: pntr in i to j => face i behind j
*/
void getAdjGraph(faces, numfaces, view, rhs, normal)
face **faces;			/* array of face pntrs, faces[0] head of lst */
int numfaces;
double *view, rhs, *normal;
{
  int numbehind, f, i, check;
  face *fpcur, *fpchk, **behind;
  extern int k_;

  if(TRUE == TRUE) {		/* if memory can be sacrificed for speed */
    /* set up huge n^2 blocked face pointer arrays for each face */
    for(f = 0; f < numfaces; f++) {
      CALLOC(faces[f]->behind, numfaces, face *, ON, AMSC);
      faces[f]->numbehind = 0;
    }      

    /* for each face, check through all faces not previously checked */
    for(fpcur = faces[0], f = 0; fpcur != NULL; fpcur = fpcur->next, f++) {
      for(fpchk = fpcur->next, numbehind = i = 0; fpchk != NULL; 
	  fpchk = fpchk->next, i++) {
	if(fpchk == fpcur) continue;	/* a face can't be behind itself */
	if((check = is1stFaceDeeper(fpcur, fpchk, view, rhs, normal))==TRUE) { 
	  fpcur->behind[(fpcur->numbehind)++] = fpchk;
	}
	else if(check == REVERSE) fpchk->behind[(fpchk->numbehind)++] = fpcur;
      }
      if(f % 20 == 0 && f != 0) {
	fprintf(stdout, "%d ", f);
	fflush(stdout);
      }
      if(f % 200 == 0 && f != 0) fprintf(stdout, "\n");
    }
  }
  else {
    CALLOC(behind, numfaces, face *, ON, AMSC);

    /* for each face, check through all the other faces */
    for(fpcur = faces[0], f = 0; fpcur != NULL; fpcur = fpcur->next, f++) {
      for(fpchk = faces[0], numbehind = i = 0; fpchk != NULL; 
	  fpchk = fpchk->next, i++) {
	if(fpchk == fpcur) continue;	/* a face can't be behind itself */
	if(is1stFaceDeeper(fpcur, fpchk, view) == TRUE) { 
	  behind[numbehind++] = fpchk;
	}
      }
      /* transfer blocked face pointers to face struct */
      fpcur->numbehind = numbehind;
      if(numbehind > 0) CALLOC(fpcur->behind, numbehind, face *, ON, AMSC);
      for(; numbehind > 0; numbehind--)
          fpcur->behind[numbehind-1] = behind[numbehind-1];
      if(f % 20 == 0 && f != 0) {
	fprintf(stdout, "%d ", f);
	fflush(stdout);
      }
      if(f % 200 == 0 && f != 0) fprintf(stdout, "\n");
    }
  }
}