zbufsort.c

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

/*
Copyright (c) 1990 Massachusetts Institute of Technology, Cambridge, MA.
All rights reserved.

This Agreement gives you, the LICENSEE, certain rights and obligations.
By using the software, you indicate that you have read, understood, and
will comply with the terms.

Permission to use, copy and modify for internal, noncommercial purposes
is hereby granted.  Any distribution of this program or any part thereof
is strictly prohibited without prior written consent of M.I.T.

Title to copyright to this software and to any associated documentation
shall at all times remain with M.I.T. and LICENSEE agrees to preserve
same.  LICENSEE agrees not to make any copies except for LICENSEE'S
internal noncommercial use, or to use separately any portion of this
software without prior written consent of M.I.T.  LICENSEE agrees to
place the appropriate copyright notice on any such copies.

Nothing in this Agreement shall be construed as conferring rights to use
in advertising, publicity or otherwise any trademark or the name of
"Massachusetts Institute of Technology" or "M.I.T."

M.I.T. MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.  By
way of example, but not limitation, M.I.T. MAKES NO REPRESENTATIONS OR
WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR
THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS OR DOCUMENTATION WILL
NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS.
M.I.T. shall not be held liable for any liability nor for any direct,
indirect or consequential damages with respect to any claim by LICENSEE
or any third party on account of or arising from this Agreement or use
of this software.
*/

#include "mulGlobal.h"
#include "zbufGlobal.h"

int cnt;			/* used in setting up the depth graph */

/*
  returns TRUE if difference of two doubles is within machine precision
*/
int diff_is_zero(num1, num2, bias)
double num1, num2, bias;
{
  double margin;

  margin = MAX(fabs(num1), fabs(num2));
  
  if(fabs(num1-num2) <= margin*MARGIN*bias + MARGIN*bias) return(TRUE);
  return(FALSE);
}
    

/*
  returns TRUE if diff. of two doubles is negative within machine precision
  - note that argument order is relevant
*/
int diff_is_negative(num1, num2, bias)
double num1, num2, bias;
{
  double margin;

  margin = MAX(fabs(num1), fabs(num2));
  
  if(num1-num2 < -margin*MARGIN*bias - MARGIN*bias) return(TRUE);
  return(FALSE);
}
    
/*
  returns the dot product
*/
double dot(vec1, vec2)
double *vec1, *vec2;
{
  int i;
  double ret = 0.0;
  
  for(i = 0; i < 3; i++) ret += vec1[i]*vec2[i];
  return(ret);
}

/*
  calculate the cross product out = in1Xin2
*/
void crossProd(out, in1, in2)
double *out, *in1, *in2;
{
  out[0] = in1[1]*in2[2] - in1[2]*in2[1];
  out[1] = in1[2]*in2[0] - in1[0]*in2[2];
  out[2] = in1[0]*in2[1] - in1[1]*in2[0];
}

/*
  returns a normal and rhs for a plane given by three points
*/
double getPlane(normal, p1, p2, p3)
double *normal, *p1, *p2, *p3;
{
  int i;
  double dot(), v12[3], v13[3], norm;

  /* set up side vectors */
  for(i = 0; i < 3; i++) {
    v12[i] = p2[i]-p1[i];
    v13[i] = p3[i]-p1[i];
  }

  /* get the normal using a cross product */
  crossProd(normal, v12, v13);
  norm = sqrt(dot(normal, normal));
  for(i = 0; i < 3; i++) normal[i] = normal[i]/norm;

  /* set up the rhs using p1 as a point on the plane */
  return(dot(normal, p1));
}

/*
  returns POS, NEG or SPLIT for which side face corners are rel to plane
*/
int whichSide(fac, facplane)
face *fac, *facplane;
{
  int i, neg, pos, zero;
  double value[MAXSIDES];      	/* holds values when subbed in plane equ */
  double dot(), margin[MAXSIDES], temp[MAXSIDES];

  /* if planes are identical, return SAME right away */
  if(fac->rhs < facplane->rhs+MARGIN && fac->rhs > facplane->rhs-MARGIN && 
     fac->normal[0] < facplane->normal[0]+MARGIN &&
     fac->normal[0] > facplane->normal[0]-MARGIN &&
     fac->normal[1] < facplane->normal[1]+MARGIN && 
     fac->normal[1] > facplane->normal[1]-MARGIN && 
     fac->normal[2] < facplane->normal[2]+MARGIN &&
     fac->normal[2] > facplane->normal[2]-MARGIN) return(SAME);

  /* sub all the points of fac into the equation for the plane of other
     - POS if result of all subs is positive or zero (normals pnt to view)
     - NEG if result of all subs is negative or zero
     - SPLIT whenever at least one positive and at least one negative
     - SAME if faces are in the same plane */
  for(i = 0; i < fac->numsides; i++) {
    temp[i] = dot(facplane->normal, fac->c[i]);
    value[i] = temp[i] - facplane->rhs;
    margin[i] = MARGIN*MAX(fabs(facplane->rhs), fabs(temp[i]));
  }
  /* count the number of each type of corner */
  for(i = neg = pos = zero = 0; i < fac->numsides; i++) {
    if(value[i] >= -margin[i] && value[i] <= margin[i]) zero++;
    else if(value[i] > 0.0) pos++;
    else if(value[i] < 0.0) neg++;
    else {
      fprintf(stderr, 
	      "\nwhichSide: strange corner, value = %g, rhs = %g, dot = %g\n",
	      value[i], facplane->rhs, temp[i]);
      exit(1);
    }
  }
  if(neg > 0 && pos == 0) return(NEG);
  else if(pos > 0 && neg == 0) return(POS);
  else if(pos > 0 && neg > 0) return(SPLIT);
  else if(zero == fac->numsides) return(SAME);
  else {
    fprintf(stderr,"\nwhichSide: face has %d corners, %d pos %d neg %d zero\n",
	    fac->numsides, pos, neg, zero);
    exit(1);
  }
}

/*
  returns TRUE if lines from1-to1 and from2-to2 intersect
  - if lines are parallel up to MARGIN, FALSE is returned
  - lines than hit at segment endpoints up to MARGIN return FALSE
*/
int doLinesIntersect(isect, from1, to1, from2, to2)
double *from1, *to1, *from2, *to2, *isect;
{
  double A[2][2], b[2];
  double det, margin, temp1, temp2, margin1, margin2;

#if DEBUGX == ON
  fprintf(stdout, "seg1 = (%g %g) (%g %g) seg2 = (%g %g) (%g %g)\n",
	  from1[0], from1[1], to1[0], to1[1], from2[0], from2[1], 
	  to2[0], to2[1]);
#endif

  /* build the 2x2 system to solve for alpha1, alpha2 in
     from1 + alpha1*(to1 - from1) = from2 + alpha2*(to2 - from2) */
  A[0][0] = to1[0]-from1[0];
  A[1][0] = to1[1]-from1[1];
  A[0][1] = from2[0]-to2[0];
  A[1][1] = from2[1]-to2[1];
  b[0] = from2[0]-from1[0];
  b[1] = from2[1]-from1[1];

  /* attempt to solve, checking carefully for cancellation in the det */
  temp1 = A[0][0]*A[1][1];
  temp2 = A[1][0]*A[0][1];
  /* margin = MARGIN*MAX(temp1, temp2); */
  det = temp1 - temp2;
  if(diff_is_zero(temp1, temp2, 1.0)) return(FALSE);	/* check for || */

  /* solve */
  /* margin1 = MARGIN*MAX(A[1][1]*b[0], A[0][1]*b[1]);
  margin2 = MARGIN*MAX(A[0][0]*b[1], A[1][0]*b[0]); */
  temp1 = A[1][1]*b[0]-A[0][1]*b[1]; /* temp = alpha*det */
  temp2 = A[0][0]*b[1]-A[1][0]*b[0];

  /* check for end hits => no crossing (just barely touching) */
  if(diff_is_zero(A[1][1]*b[0], A[0][1]*b[1], 1.0) ||
     diff_is_zero(A[1][1]*b[0], A[0][1]*b[1]+det, 1.0)) return(FALSE);
  if(diff_is_zero(A[0][0]*b[1], A[1][0]*b[0], 1.0) ||
     diff_is_zero(A[0][0]*b[1], A[1][0]*b[0]+det, 1.0)) return(FALSE);

  /* set up intersection point using isect = from1 + alpha1*(to1 - from1) */
  isect[0] = from1[0] + temp1*(to1[0] - from1[0])/det;
  isect[1] = from1[1] + temp1*(to1[1] - from1[1])/det;
  isect[2] = 0.0;		/* not really necessary */

  /* check for alphas between 0 and 1 (temps btwn 0 and det) */
  if(det > 0.0) {
    if(0.0 < temp1 && temp1 < det && 0.0 < temp2 && temp2 < det) return(TRUE);
    else return(FALSE);
  }
  else {
    if(det < temp1 && temp1 < 0.0 && det < temp2 && temp2 < 0.0) return(TRUE);
    else return(FALSE);
  }
}
  
/*
  returns TRUE if one panel is completely inside the other
  - does careful checks when vertices on one face are on verts, sides of other
*/
int face_is_inside(corners1, ccnt1, corners2, ccnt2, com_pnt)
double **corners1, **corners2, *com_pnt;
int ccnt1, ccnt2;
{
  int i, j, k, n, ccnt, zeros, pos, neg, ncnt;
  double *sfrom, *sto, margin1, margin2, **refcor, **curcor;
  double innerpnt[2], side[2], pnt[2];
  double temp;

  /* do cross products between sides of one face and vectors to first point
     of the other---must get same sign or zero for all if one inside other 
     - assumes faces have no sides that intersect 
     - do assuming both sides are outside to be sure to catch all cases */
  for(refcor = corners1, curcor = corners2, ccnt = ccnt1, ncnt = ccnt2, n = 0;
      n < 2; 
      refcor = corners2, curcor = corners1, ccnt = ccnt2, ncnt = ccnt1, n++) {
    zeros = pos = neg = 0;	/* cross-product value counts */
    for(i = 0; i < ccnt; i++) {	/* loop on "outer" face sides */
      if(i == ccnt-1) j = 0;
      else j = i+1;
      /* figure the center point of the non-traversed panel */
      innerpnt[0] = innerpnt[1] = 0.0;
      for(k = 0; k < ncnt; k++) {
	innerpnt[0] += curcor[k][0];
	innerpnt[1] += curcor[k][1];
      }
      innerpnt[0] /= (double)ncnt;
      innerpnt[1] /= (double)ncnt;
      /* figure side and point vector */
      side[0] = refcor[j][0] - refcor[i][0];
      side[1] = refcor[j][1] - refcor[i][1];
      pnt[0] = innerpnt[0] - refcor[i][0];
      pnt[1] = innerpnt[1] - refcor[i][1];
      margin1 = MARGIN*MAX(fabs(innerpnt[0]), fabs(refcor[i][0]));
      margin2 = MARGIN*MAX(fabs(innerpnt[1]), fabs(refcor[i][1]));
      /* figure cross-product---if a cross-product would be zero, skip it */
      if(diff_is_zero(innerpnt[0], refcor[i][0], 1.0)
	 && diff_is_zero(innerpnt[1], refcor[i][1], 1.0)) {
	zeros++;
	continue;
      }
      /* now check if vertex of "inside" face is on a side, not on corner */
      if(diff_is_zero(side[0]*pnt[1], pnt[0]*side[1], 1.0)) { /* is X-prod 0 */
	zeros++;		/*       due to vertex on side (not corner) */
      }
      else if(diff_is_negative(side[0]*pnt[1], pnt[0]*side[1], 1.0)) {
	neg++;
      }
      else pos++;		/* solid positive cross-product */
    }
    /* see if negative and positive are not mixed together => full overlap */
    if(((neg == 0 && pos != 0) || (neg != 0 && pos == 0)) && zeros == 0) {
      for(k = 0; k < 3; k++) com_pnt[k] = innerpnt[k];
      return(TRUE);
    }
  }

  /* wasn't any full overlap in either case => no overlap */
  return(FALSE);
	
}
  
#if 1 == 0			/* old is1stFaceDeeper */
/*
  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)
face *fac, *facref;
double *view;
{
  int i, j, k, olap[2], is_overlap;
  static double **cproj = NULL;	/* corners of fac in facref's plane */
  static double **cref = NULL;	/* 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, temp1;
  double margin1;
  double *cfr, *ctr, *cff, *ctf, avg[3];
  int all_pos, all_neg, intersect;

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

  /* allocate for local arrays on first call */
  if(cproj == NULL) {
    CALLOC(cproj, MAXSIDES, double *, ON, AMSC);
    CALLOC(cref, MAXSIDES, double *, ON, AMSC);
    for(i = 0; i < MAXSIDES; i++) {
      CALLOC(cproj[i], 3, double, ON, AMSC);
      CALLOC(cref[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 */
  temp = temp1 = margin = 0.0;
  for(i = 0; i < 3; i++) {
    temp1 = fac->normal[i]-facref->normal[i];
    margin = MAX(margin, fabs(fac->normal[i]));
    margin = MAX(margin, fabs(facref->normal[i]));
    temp += (temp1*temp1);
  }
  temp = sqrt(temp);		/* get norm of difference of normals */
  margin *= MARGIN;		/* get norm of difference margin */
  /* check rhs and normal equivalence (panels in same plane) */
  if(diff_is_zero(fac->rhs, facref->rhs, 1.0) 
     && -margin <= temp && temp <= margin) {
    return(FALSE);
  }

  /* 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;
    /* test fails if v-c is perpendicular to n */
    if(temp > -margin && temp < margin) {
      fprintf(stderr, 
	      "is1stFaceDeeper: Warning, view-corner in view plane (?)\n");
      return(FALSE);
    }
    /* get alpha as in c + alpha*(v-c) = c' */
    alpha[i] = (facref->rhs - dot(facref->normal, fac->c[i]))/temp;
    for(j = 0; j < 3; j++)	/* get c' */
	cproj[i][j] = (1.0-alpha[i])*fac->c[i][j]+alpha[i]*view[j];
  }