zbufproj.c

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

      miny = MAX(lines[i]->from[1], lines[i]->to[1]);
    }
    else {
      minx = MIN(minx, lines[i]->from[0]);
      minx = MIN(minx, lines[i]->to[0]);
      miny = MIN(miny, lines[i]->from[1]);
      miny = MIN(miny, lines[i]->to[1]);
    }
  }
  /* include axes if they will be printed */
  for(i = 0; i < 7 && x_ == TRUE; i++) {
    for(j = 0; j < 2; j++) {
      minx = MIN(minx, axes[i][j][0]);
      miny = MIN(miny, axes[i][j][1]);
    }
  }
    
  /* transformation makes minx = offset[0], miny = offset[1] */
  trans[0] = offset[0]-minx;
  trans[1] = offset[1]-miny;
  for(i = 0; i < numfaces; i++) {
    for(j = 0; j < faces[i]->numsides; j++) {
      (faces[i]->c)[j][0] += trans[0];
      (faces[i]->c)[j][1] += trans[1];  
    }  
  }
  /* xform lines */
  for(i = 0; i < numlines; i++) {
    (lines[i]->from)[0] += trans[0];
    (lines[i]->from)[1] += trans[1];  
    (lines[i]->to)[0] += trans[0];
    (lines[i]->to)[1] += trans[1];  
  }
  for(i = 0; i < 7 && x_ == TRUE; i++) {  /* include axes if to be printed */
    for(j = 0; j < 2; j++) {
      axes[i][j][0] += trans[0];
      axes[i][j][1] += trans[1];  
    }  
  }
}

/*
  scales the 2d image so that it's offsetted and scale*IMAGEX x scale*IMAGEY
  - the global defines OFFSETX/Y set the offset
  - assumes smallest x and y coordinates are zero (called after makePos())
*/
void scale2d(faces, numfaces, lines, numlines, scale, offset)
face **faces;
line **lines;
int numfaces, numlines;
double scale, *offset;
{
  int i, j;
  double ymax, xmax, xscale, yscale, master;
  extern double ***axes;
  extern int x_;

  /* find ymax, xmax */
  for(i = 0, xmax = ymax = 0.0; i < numfaces; i++) {
    for(j = 0; j < faces[i]->numsides; j++) {
      xmax = MAX(xmax, faces[i]->c[j][0]);
      ymax = MAX(ymax, faces[i]->c[j][1]);
    }
  }
  /* include lines */
  for(i = 0; i < numlines; i++) {
    xmax = MAX(xmax, lines[i]->to[0]);
    ymax = MAX(ymax, lines[i]->to[1]);
    xmax = MAX(xmax, lines[i]->from[0]);
    ymax = MAX(ymax, lines[i]->from[1]);
  }
  for(i = 0; i < 7 && x_ == TRUE; i++) { /* if axes to be prnted */
    for(j = 0; j < 2; j++) {
      xmax = MAX(xmax, axes[i][j][0]);
      ymax = MAX(ymax, axes[i][j][1]);
    }
  }

  if(xmax <= MARGIN || ymax <= MARGIN) {
    fprintf(stderr, "\nscale2d: strange xmax = %g or ymax = %g\n",
	    xmax, ymax);
    exit(1);
  }

  /* find the x and y scales that would make those dimensions dead on */
  xscale = scale*IMAGEX/xmax;
  yscale = scale*IMAGEY/ymax;
  
  /* figure the final scaling and apply it - do the offsetting too */
  master = MIN(xscale, yscale);
  for(i = 0; i < numfaces; i++) {
    for(j = 0; j < faces[i]->numsides; j++) {
      faces[i]->c[j][0] *= master; 
      faces[i]->c[j][1] *= master;
      faces[i]->c[j][0] += OFFSETX; 
      faces[i]->c[j][1] += OFFSETY;
    }
  }
  for(i = 0; i < numlines; i++) {
    lines[i]->to[0] *= master; 
    lines[i]->to[1] *= master;
    lines[i]->to[0] += OFFSETX; 
    lines[i]->to[1] += OFFSETY;
    lines[i]->from[0] *= master; 
    lines[i]->from[1] *= master;
    lines[i]->from[0] += OFFSETX; 
    lines[i]->from[1] += OFFSETY;
  }
  for(i = 0; i < 7 && x_ == TRUE; i++) { /* do axes if needed */
    for(j = 0; j < 2; j++) {
      axes[i][j][0] *= master; 
      axes[i][j][1] *= master;
      axes[i][j][0] += OFFSETX; 
      axes[i][j][1] += OFFSETY;
    }
  }
}

/*
  returns the center of a rectangular prism contianing all the face corners
*/
double *getAvg(faces, numfaces, lines, numlines, flag)
face **faces;
line **lines;
int numfaces, numlines;
int flag;			/* ON => include axes */
{
  double *avg, max[3], min[3];
  int i, j, k;
  extern double ***axes;
  extern int x_;

  CALLOC(avg, 3, double, ON, AMSC);

  /* check the faces' coordinates */
  for(i = 0; i < numfaces; i++) {
    for(j = 0; j < faces[i]->numsides; j++) {
      if(i == 0 && j == 0) {
	for(k = 0; k < 3; k++) max[k] = min[k] = (faces[0]->c)[0][k];
      }
      else {
	for(k = 0; k < 3; k++) {
	  max[k] = MAX(max[k], faces[i]->c[j][k]);
	  min[k] = MIN(min[k], faces[i]->c[j][k]);
	}
      }
    }
  }

  /* check the lines' coordinates */
  for(i = 0; i < numlines; i++) {
    if(i == 0 && numfaces == 0) {
      for(k = 0; k < 3; k++) {
	max[k] = MAX((lines[0]->from)[k], (lines[0]->to)[k]);
	min[k] = MIN((lines[0]->from)[k], (lines[0]->to)[k]);
      }
    }
    else {
      for(k = 0; k < 3; k++) {
	max[k] = MAX(max[k], lines[0]->from[k]);
	max[k] = MAX(max[k], lines[0]->to[k]);
	min[k] = MIN(min[k], lines[0]->from[k]);
	min[k] = MIN(min[k], lines[0]->to[k]);
      }
    }
  }

  /* if axes are specified, they must be included */
  for(i = 0; i < 7 && x_ == TRUE && flag == ON; i++) {
    for(j = 0; j < 2; j++) {
      for(k = 0; k < 3; k++) {
	max[k] = MAX(max[k], axes[i][j][k]);
	min[k] = MIN(min[k], axes[i][j][k]);
      }
    }
  }

  /* average each coordinates min an max to get center of picture */
  for(k = 0; k < 3; k++) avg[k] = (max[k] + min[k])/2;

  return(avg);
}

/*
  returns the radius of the smallest sphere with center at avg that
  encloses all the faces given in faces, all lines and axes, if required
*/
double getSphere(avg, faces, numfaces, lines, numlines)
face **faces;
line **lines;
double *avg;
int numfaces, numlines;
{
  double radius = 0.0, dot(), temp[3];
  int i, j, k, l;
  extern int x_;
  extern double ***axes;

  /* loop through all the points, save the farthest away */
  for(i = 0; i <  numfaces; i++) {
    for(j = 0; j < faces[i]->numsides; j++) {
      for(k = 0; k < 3; k++) temp[k] = avg[k] - faces[i]->c[j][k];
      radius = MAX(radius, dot(temp, temp));
    }
  }

  /* loop through all the line points, save the farthest away */
  for(i = 0; i < numlines; i++) {
    for(k = 0; k < 3; k++) temp[k] = avg[k] - lines[i]->to[k];
    radius = MAX(radius, dot(temp, temp));
    for(k = 0; k < 3; k++) temp[k] = avg[k] - lines[i]->from[k];
    radius = MAX(radius, dot(temp, temp));
  }

  /* if axes are specified, they must be included */
  for(i = 0; i < 7 && x_ == TRUE; i++) {
    for(j = 0; j < 2; j++) {
      for(k = 0; k < 3; k++) temp[k] = avg[k] - axes[i][j][k];
      radius = MAX(radius, dot(temp, temp));
    }
  }

  return(sqrt(radius));
}

/* 
  get the normal to the image plane, adjust view point to be 2*radius away
  from object center - view coordinates are azimuth and elevation
  relative to object center at first
  - also change axes lengths so that they are smaller than view distance
    (ie so they won't cross the view plane)
*/
double getNormal(normal, radius, avg, view, distance)
double *normal, radius, *avg, *view, distance;
{
  int i, k, j, axes_too_big, first;
  double rhs, dot(), norm, anorm, max_anorm;
  extern int x_;
  extern double ***axes;

  /* figure the normal to the view plane 
     - remember view is rel to avg
     - use view[0] as azimuth relative to poitive x direction
     - use view[1] as elevation relative to positive z direction */
  normal[0] = normal[1] = sin(M_PI*view[1]/180.0);
  normal[0] *= cos(M_PI*view[0]/180.0);	/* x = sin(beta)cos(alpha) */
  normal[1] *= sin(M_PI*view[0]/180.0);	/* y = sin(beta)sin(alpha) */
  normal[2] = cos(M_PI*view[1]/180.0); /* z = cos(beta) */

  /* figure the rhs of the plane's equation (use view for temp storage) */
  for(i = 0; i < 3; i++) view[i] = (1.0+distance/2.0)*radius*normal[i]+avg[i];
  rhs = dot(normal, view);

  /* adjust view point - change to absolute coordinates */
  for(i = 0; i < 3; i++) view[i] = (1.0+distance)*radius*normal[i] + avg[i];

  if(x_ == FALSE) {
    /* adjust the axes if needed to keep from going beyond the view plane */
    /*    get norm of view-avg */
    for(norm = 0.0, i = 0; i < 3; i++) 
	norm += (view[i]-avg[i])*(view[i]-avg[i]);
    /*    check sizes of current axes */
    axes_too_big = FALSE;
    first = TRUE;
    for(i = 0; i < 7; i++) {
      anorm = 0.0;
      for(k = 0; k < 3; k++) {
	anorm += (axes[i][0][k]-axes[i][1][k])*(axes[i][0][k]-axes[i][1][k]);
      }
      if(anorm >= norm) {
	axes_too_big = TRUE;
	if(first) {
	  first = FALSE;
	  max_anorm = anorm;
	}
	else max_anorm = MAX(max_anorm, anorm);
      }
    }
    /*    adjust if too big */
    if(axes_too_big) {
      /* make the biggest anorm = half the view-avg norm */
      max_anorm = sqrt(max_anorm); /* the biggest norm */
      anorm = sqrt(norm)/max_anorm; /* the scale factor */
      anorm /= 2.0;
      for(i = 0; i < 7; i++) {
	for(j = 0; j < 2; j++) {
	  for(k = 0; k < 3; k++) {
	    axes[i][j][k] *= anorm;
	  }
	}
      }
    }
  }
  return(rhs);
}