zbufsort.c

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

  /* 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 */
  }

  /* up to here is identical to isThereBoxOverlap() */
  /* for each side of facref, see if there is an intersect. w/ all fac sides */
#if DEBUGX == ON
  fprintf(stdout, "Is face %d behind face %d?\n", fac->index, facref->index);
#endif
  is_overlap = FALSE;
  for(i = 0; i < facref->numsides && !is_overlap; i++) {
    cfr = cref[i];
    if(i == facref->numsides - 1) ctr = cref[0];
    else ctr = cref[i+1];
    for(j = 0; j < fac->numsides && !is_overlap; j++) {
      cff = cproj[j];
      if(j == fac->numsides - 1) ctf = cproj[0];
      else ctf = cproj[j+1];
      if((intersect = doLinesIntersect(cfr, ctr, cff, ctf)) == TRUE)
	  is_overlap = TRUE;
#if DEBUGX == ON
      fprintf(stdout, "doLinesIntersect returned %d\n", intersect);
#endif
    }
  }
#if XOVTST == ON		/* do either this or face_is_inside() below */
  /* check for overlap with lines across face */
  if(facref->numsides == 4) {
    for(i = 0; i < 2; i++) {
      if(i == 0) {
	cfr = cref[0];
	ctr = cref[2];
      }
      else {
	cfr = cref[1];
	ctr = cref[3];
      }
      for(j = 0; j < fac->numsides; j++) {
	cff = cproj[j];
	if(j == fac->numsides - 1) ctf = cproj[0];
	else ctf = cproj[j+1];
	if((intersect = doLinesIntersect(cfr, ctr, cff, ctf)) == TRUE)
	    is_overlap = TRUE;
#if DEBUGX == ON
	fprintf(stdout, "doLinesIntersect returned %d\n", intersect);
#endif
      }
    }
  }
  else if(facref->numsides == 3) {
    avg[0] = (cref[0][0] + cref[1][0] + cref[2][0])/3.0;
    avg[1] = (cref[0][1] + cref[1][1] + cref[2][1])/3.0;
    avg[2] = (cref[0][2] + cref[1][2] + cref[2][2])/3.0;
    cfr = avg;
    for(i = 0; i < 3; i++) {
      ctr = cref[i];
      for(j = 0; j < fac->numsides; j++) {
	cff = cproj[j];
	if(j == fac->numsides - 1) ctf = cproj[0];
	else ctf = cproj[j+1];
	if((intersect = doLinesIntersect(cfr, ctr, cff, ctf)) == TRUE)
	    is_overlap = TRUE;
#if DEBUGX == ON
	fprintf(stdout, "doLinesIntersect returned %d\n", intersect);
#endif
      }
    }
  }
  else {
    fprintf(stderr, "isThereProjOverlap: can't handle %d side panel\n",
	    facref->numsides);
    exit(1);
  }
#else
  /* no sides overlap---check if one face completely obscures the other */
  if(!is_overlap) {
    is_overlap = face_is_inside(cref, facref->numsides, cproj, fac->numsides);
  }
#endif				/* XOVTST == ON */

  /* if no overlap, no edge in graph in any case */
  if(!is_overlap) return(FALSE);

  /* return TRUE only if fac is deeper than facref */
  /* return REVERSE only if facref is deeper than fac */
  /* if all alphas are positive or zero, fac is deeper (note use of margin) */
  /* if all alphas are negative or zero, facref is deeper */
  /* otherwise the test is inconclusive */
  all_pos = all_neg = TRUE;
  for(i = 0; i < fac->numsides; i++) {
    if(alpha[i] < -margin) all_pos = FALSE;
    else if(alpha[i] > margin) all_neg = FALSE;
  }

  if(all_pos && all_neg) {
    fprintf(stderr, 
	 "\nis1stFaceDeeper: Warning, all projection coefficients are zero\n");
    return(FALSE);
  }
  else if(all_pos) return(TRUE);
  else if(all_neg) return(REVERSE);
  else return(FALSE);
}

#else				/* new is1stFaceDeeper--proj to view plane */

/*
  returns TRUE if fac is deeper than facref and facref overlaps fac
  returns FALSE if facref has no overlap with fac
  returns REVERSE if facref is deeper than fac and fac overlaps facref
  - checks for intersections between each line of fac and all sides of facref
  - also checks for complete overlap (one face inside the other)
*/
int is1stFaceDeeper(fac, facref, view, rhs, normal)
face *fac, *facref;
double *view, rhs, *normal;
{
  int i, j, k, olap[2], is_overlap, isect_cnt;
  static double ***cproj = NULL;	/* corners of faces in view plane */
  double alpha[2][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, temp1;
  double margin1;
  double *cfr, *ctr, *cff, *ctf, avg[3], origin[3];
  double isect_avg[3], isect[3]; /* intersection points */
  double alpha_fac, alpha_facref, alpha_origin, alpha_x;
  int all_pos, all_neg, intersect, same_normal;
  face *curf;

  if(fac->index == 8 && facref->index == 20
     || fac->index == 20 && facref->index == 8) {
    i = i + 0;
  }

  /* allocate for local arrays on first call */
  if(cproj == NULL) {
    CALLOC(cproj, 2, double **, ON, AMSC);
    for(k = 0; k < 2; k++) {
      CALLOC(cproj[k], MAXSIDES, double *, ON, AMSC);
      for(i = 0; i < MAXSIDES; i++) {
	CALLOC(cproj[k][i], 3, double, ON, AMSC);
      }
    }
  }

  /* figure if panels are in the same plane
     - if they are, they cant overlap if this is a legal discretization
       so return false */
#if 1 == 0
  temp = temp1 = margin = 0.0;
  for(i = 0; i < 3; i++) {
    temp1 = fac->normal[i]-facref->normal[i];
    temp1 = temp1*temp1;
    margin = MAX(margin, temp1);
    temp += temp1;
  }
  temp = sqrt(temp);		/* get norm of difference of normals */
  margin *= (MARGIN*PARMGN);	/* get norm of difference margin */
#endif
  same_normal = TRUE;
  for(i = 0; i < 3 && same_normal; i++) {
    if(!diff_is_zero(fac->normal[i], facref->normal[i], PARMGN)) 
	same_normal = FALSE;
  }
    
  /* check rhs and normal equivalence (panels in same plane) */
  if(diff_is_zero(fac->rhs, facref->rhs, PARMGN) && same_normal) return(FALSE);

  /* find projections of fac and facref corners onto view plane rel to view */
  for(curf = fac, k = 0; k < 2; k++, curf = facref) {
    for(i = 0; i < curf->numsides; i++) {
      for(j = 0; j < 3; j++) tvec[j] = view[j] - curf->c[i][j]; /* get v-c */
      temp = dot(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[k][i] = (rhs - dot(normal, curf->c[i]))/temp;
      for(j = 0; j < 3; j++)	/* get c' */
     /*cproj[k][i][j] = (1.0-alpha[k][i])*curf->c[i][j]+alpha[k][i]*view[j];*/
	  cproj[k][i][j] = curf->c[i][j] + alpha[k][i]*tvec[j];
    }
  }

#if 1 == 0
  /* project origin and (1 0 0) into view plane */
  /* live dangerously---no perpendicular check */
  temp = dot(normal, view);
  if(temp == 0.0) {
    fprintf(stderr, 
	"is1stFaceDeeper: origin-view vector perpendicular to view plane\n");
  }
  /* get projection alpha */
  alpha_origin = rhs/temp;
  for(j = 0; j < 3; j++) origin[j] = alpha_origin*view[j];

  /* project (1 0 0) */
  for(j = 0; j < 3; j++) tvec[j] = view[j];
  tvec[0] -= 1.0;
  if((temp = dot(normal, tvec)) == 0.0) {
    fprintf(stderr,
	    "is1stFaceDeeper: x-view vector perpendicular to view plane\n");
  }
  /* get projection alpha */
  tvec1[0] = 1.0; tvec1[1] = tvec1[2] = 0.0;
  alpha_x = (rhs - dot(normal, tvec1))/temp;
  /* get x-axis projection */
  for(j = 0; j < 3; j++) x[j] = tvec1[j] + alpha_x*tvec[j] - origin[j];
  temp = sqrt(dot(x, x));
  for(j = 0; j < 3; j++) x[j] /= temp;
  /* y = zXx */
  crossProd(y, normal, x);
#endif

#if 1 == 1
  /* figure x and y coordinates in view plane (normal is z coordinate) */
  /* x = c0-c1 projections from fac proj. always (should never be 0 len) */
  for(j = 0; j < 3; j++) origin[j] = cproj[0][0][j]; /* to stop overwrites */
  for(j = 0; j < 3; j++) x[j] = cproj[0][1][j] - origin[j];
  temp = sqrt(dot(x, x));
  for(j = 0; j < 3; j++) x[j] /= temp; /* normalize */
  /* y = zXx */
  crossProd(y, normal, x);
#endif
  
  /* project all face corner projections onto new x and y coordinates
     - cproj[0][0] plays the role of origin */
  for(curf = fac, k = 0; k < 2; k++, curf = facref) {	/* loop on faces */
    for(i = 0; i < curf->numsides; i++) {
      for(j = 0; j < 3; j++) tvec1[j] = cproj[k][i][j] - origin[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 = rhs */
      for(j = 0; j < 3; j++) cproj[k][i][j] = tvec[j]; /* xfer */
    }
  }
  for(j = 0; j < 3; j++) cproj[0][0][j] = 0.0; /* set origin explicitly zero */

  /* for each side of facref, see if there is an intersect. w/ all fac sides */
#if DEBUGX == ON
  fprintf(stdout, "Is face %d behind face %d?\n", fac->index, facref->index);
#endif
  is_overlap = FALSE;
  isect_cnt = 0;
  isect_avg[0] = isect_avg[1] = isect_avg[2] = 0.0;
  for(i = 0; i < facref->numsides; i++) {
    cfr = cproj[1][i];
    if(i == facref->numsides - 1) ctr = cproj[1][0];
    else ctr = cproj[1][i+1];
    for(j = 0; j < fac->numsides; j++) {
      cff = cproj[0][j];
      if(j == fac->numsides - 1) ctf = cproj[0][0];
      else ctf = cproj[0][j+1];
      if((intersect = doLinesIntersect(isect, cfr, ctr, cff, ctf)) == TRUE) {
	isect_cnt++;
	for(k = 0; k < 3; k++) isect_avg[k] += isect[k];
	is_overlap = TRUE;
      }
#if DEBUGX == ON
      fprintf(stdout, "doLinesIntersect returned %d\n", intersect);
#endif
    }
  }
#if XOVTST == ON		/* do either this or face_is_inside() below */
  /* check for overlap with lines across face */
  if(facref->numsides == 4) {
    for(i = 0; i < 2; i++) {
      if(i == 0) {
	cfr = cproj[1][0];
	ctr = cproj[1][2];
      }
      else {
	cfr = cproj[1][1];
	ctr = cproj[1][3];
      }
      for(j = 0; j < fac->numsides; j++) {
	cff = cproj[0][j];
	if(j == fac->numsides - 1) ctf = cproj[0][0];
	else ctf = cproj[0][j+1];
	if((intersect = doLinesIntersect(cfr, ctr, cff, ctf)) == TRUE)
	    is_overlap = TRUE;
#if DEBUGX == ON
	fprintf(stdout, "doLinesIntersect returned %d\n", intersect);
#endif
      }
    }
  }
  else if(facref->numsides == 3) {
    avg[0] = (cproj[1][0][0] + cproj[1][1][0] + cproj[1][2][0])/3.0;
    avg[1] = (cproj[1][0][1] + cproj[1][1][1] + cproj[1][2][1])/3.0;
    avg[2] = (cproj[1][0][2] + cproj[1][1][2] + cproj[1][2][2])/3.0;
    cfr = avg;
    for(i = 0; i < 3; i++) {
      ctr = cproj[1][i];
      for(j = 0; j < fac->numsides; j++) {
	cff = cproj[0][j];
	if(j == fac->numsides - 1) ctf = cproj[0][0];
	else ctf = cproj[0][j+1];
	if((intersect = doLinesIntersect(isect, cfr, ctr, cff, ctf)) == TRUE)
	    is_overlap = TRUE;
#if DEBUGX == ON
	fprintf(stdout, "doLinesIntersect returned %d\n", intersect);
#endif
      }
    }
  }
  else {
    fprintf(stderr, "isThereProjOverlap: can't handle %d side panel\n",
	    facref->numsides);
    exit(1);
  }
#else
  /* no sides overlap---check if one face completely obscures the other */
  if(!is_overlap) {
    is_overlap 
	= face_is_inside(cproj[1], facref->numsides, cproj[0], fac->numsides,
			 isect_avg);
  }
  else {
    /* figure average intersect point---should be inside overlap region */
    for(k = 0; k < 3; k++) isect_avg[k] /= (double)isect_cnt;
    if(isect_cnt % 2 != 0) {
      k = k + 0;
    }
  }
#endif				/* XOVTST == ON */

  /* if no overlap, no edge in graph in any case */
  if(!is_overlap) return(FALSE);

  /* return TRUE only if fac is deeper than facref */
  /* return REVERSE only if facref is deeper than fac */
  /* project average point back along view direction to fac, facref planes */
  /*   first convert average point back to 3D */
  isect[0] = isect_avg[0]*x[0] + isect_avg[1]*y[0];
  isect[1] = isect_avg[0]*x[1] + isect_avg[1]*y[1];
  isect[2] = isect_avg[0]*x[2] + isect_avg[1]*y[2];
  /*   add in origin to get point in absolute coordinates */
  for(k = 0; k < 3; k++) isect[k] += origin[k];
  /*   figure projection back to fac */
  for(k = 0; k < 3; k++) isect[k] -= view[k];
  alpha_fac = (fac->rhs - dot(fac->normal, view))/dot(fac->normal, isect);
  alpha_facref 
      = (facref->rhs - dot(facref->normal, view))/dot(facref->normal, isect);