wzpoint.cxx

Delaunay三角形的网格剖分程序

#include <stdio.h>
#include "wzpoint.hxx"
#include "wzoutput.hxx"

wzIndex	wzPoint::Dim		= wzPointDim;
wzIndex	wzPoint::DimFull  	= wzPointDim; // wzPointDimFull;
wzIndex	wzPoint::DimFloat 	= 0; // wzPointDimFull-wzPointDim;
wzIndex	wzPoint::DimIndex 	= wzPointDimIndex;
wzFloat	wzPoint::Epsilon     	= 1.e-6;
//const wzIndex	wzPoint::XOffset  	= wzOffsetOf(wzPoint,X);
const wzIndex  wzInvalidPointMarker	= 0xffff;

int wzAddPoints(wzPointList& list, FILE* file, int dim, int lines);
int wzAddPoints(wzPointList& p, char* f);

// This function opens the file and assumes that it contains a point
// list of the following type:
// -  may be some header lines not beginning with real numbers
// -  lines of real numbers, for each point one line, and all coordinates
//    of the point given on the same line. 
//    All lines are assumed to contain the same number of coordinates.
// -  may be some additional data with different format.
// The header lines will be ommitted, until a line is found which looks
// like a correct point description line. The dimension will be
// derived from the number of reals in this line.
// A problem are header lines beginning with integers. 
// They may be misinterpreted as a point description line.

wzFloat  wzPoint::vector2(const wzPoint& p0, const wzPoint& p1) const
{
  wzFloat       x0  =  (p0[0]-X[0]), y0= (p1[0]-X[0]); 
  wzFloat       x1  =  (p0[1]-X[1]), y1= (p1[1]-X[1]); 
  wzFloat       x2  =  (p0[2]-X[2]), y2= (p1[2]-X[2]);
  wzFloat 	z0  =  x1*y2 - x2*y1;
  wzFloat 	z1  =  x2*y0 - x0*y2;
  wzFloat 	z2  =  x0*y1 - x1*y0;
  return	z0*z0+z1*z1+z2*z2;
}

void wzPoint::affine(wzFloat p, const wzPoint& p0, const wzPoint& p1)
{register wzFloat q=1-p;
 for(int i=0;i<DimFull;i++){
        X[i] = p*p0[i] + q*p1[i];
 }
}

wzSegment      wzNoSegment	 (0,0);
wzSegment      wzRegionOutside (1,0);
wzSegment      wzRegionDefault (2,0);
wzSegment      wzFaceDefault   (1,1);
wzSegment      wzEdgeDefault   (1,2);
wzSegment      wzNodeDefault   (1,3);

void wzPoint::print() const
{
  wzOutput& pout = wzOutput::Default;
  if(inRegion())	{pout("region <> [<10>,<10>,<10>]\n"),(wzIndex)region();}
  else if(onFace())	{pout("face <> [<10>,<10>,<10>]\n"),(wzIndex)face();}
  else if(onEdge())	{pout("edge <> [<10>,<10>,<10>]\n"),(wzIndex)edge();}
  else if(isVertex())	{pout("vertex <> [<10>,<10>,<10>]\n"),(wzIndex)vertex();}
  pout<<x(),y(),z();
}

/*
wzPointList::wzPointList()
:points(sizeof(wzPoint),wzInvalidPointMarker)
,point(points.base)
,spaceDim(0),fullDim(0)
{;}
*/

wzPointList::wzPointList(char *filename)
:wzRange(sizeof(wzPoint),wzInvalidPointMarker)
,Point(base)
,Near(*(wzRange*)this,sizeof(wzIndex))
,spaceDim(0)
{wzAddPoints(*this,filename);}

int wzAddPoints(wzPointList& p, FILE* f, int dim, int fulldim, int lines=0);

// it is assumed that the file f is already open for reading.
// The header is assumed to be read, the information from the
// header is assumed to be localized in dim and lines.
// dim describes the number of reals in each line,
// lines the number of lines, lines=0 means up to the end of the file.
// Returns the number of successfully added points.

int wzAddPoints(wzPointList& p, char *filename)
{int rc,rd,rn,i,headlines,notthefirst,d[6];
 float f[6];
 const wzIndex buflen=256;
 char buffer[buflen];
 FILE *file = fopen(filename,"r"); 
 if(!file){throw wzFileOpenError();}
 headlines = 0; notthefirst=0;
 while(fgets(buffer,buflen,file)){
	rc = sscanf(buffer,"%f %f %f %f %f %f",
			&f[0],&f[1],&f[2],&f[3],&f[4],&f[5]);
	if(rc==0){
		headlines++; continue;
	}
// now we have to guess if this may nonetheless a header line.
tests:
// Let's see if they are all integers. 
// Non-integer data are not typical for a header line.
	rd = sscanf(buffer,"%d %d %d %d %d %d",
			&d[0],&d[1],&d[2],&d[3],&d[4],&d[5]);
	if(rd < rc) goto found;
// let's compare the number from this line with the next line:
	fgets(buffer,buflen,file);
	rn = sscanf(buffer,"%f %f %f %f %f %f",
			&f[0],&f[1],&f[2],&f[3],&f[4],&f[5]);
	if(rn != rc){
		headlines++; rc=rn; notthefirst = 1; goto tests;
	}
	rd = sscanf(buffer,"%d %d %d %d %d %d",
			&d[0],&d[1],&d[2],&d[3],&d[4],&d[5]);
	if(rd < rc){
// the next line is the point line. But was the previous? 
// The same number of components, but all integer. 

// If there has been already a header line with numbers before, probably yes.
		if(notthefirst) goto found;
// else, if the number of components is greater 1, probably yes.
		if(rc>1) goto found;
// else probably not.
		headlines++; goto found;
	}
// now, the next line has the same number of components, and all are integer.
// if the number of components is 1, probably it is a header line.
	if(rc==1){
		headlines++;  rc=rn; notthefirst = 1; goto tests;
	}
	goto found;
 }
found:
 fclose(file);
 file = fopen(filename,"r");
 for(i=0;i<headlines;i++) fgets(buffer,buflen,file);
 rc = wzAddPoints(p,file,rc,rc,0);
 fclose(file);
 return rc;
}

int wzAddPoints(wzPointList& list, FILE* file, int dim, int fulldim, int lines) 
{
 const int buflen=256;
 char buffer[buflen];
 float x[6];
 int p,i,rc;
 wzPoint *P;
 int ll = lines;
 wzAssert(fulldim>=dim);
 if(list.fullDim){
   if(list.fullDim!=fulldim) goto readerror;
 }else{
   list.fullDim = fulldim;
 }
 if(list.spaceDim){
   if(list.spaceDim!=dim) goto readerror;
 }else{
   list.spaceDim = dim;
 }
 while(fgets(buffer,buflen,file)){
    	rc = sscanf(buffer,"%f %f %f %f %f %f",
				&x[0],&x[1],&x[2],&x[3],&x[4],&x[5]);
	if(rc!=fulldim){
	  if(ll<0) break;
	  goto readerror;
	}
	p = list.create();
	P = &list.Point[p];
	list.Near[p] = 0;
	for(i=0;i<fulldim;i++){
		(*P)[i] = x[i];
	}
	if((--ll)==0) break;
 }
 if(list.spaceDim>3) goto readerror; 
 return wzSuccess;
readerror:
 throw wzFileReadError();
}

void wzPointList::setFloatValues(wzIndex fdim)
{
  fDim = fdim; fullDim = spaceDim+fDim;
  wzPoint p; p.setFloatValues(fDim);
}


void wzPointList::computeBox()
{
  int i,j;
  wzPoint* P;
  for(j=0;j<spaceDim;j++){
    min[j] = wzInfty;
    max[j] = -wzInfty;
  }
  for(i=1;i<=last();i++){
    P = &(Point[i]);
    for(j=0;j<spaceDim;j++){
      if((*P)[j]<min[j])  min[j] = (*P)[j];
      if((*P)[j]>max[j])  max[j] = (*P)[j];
    }
  }
  for(j=0;j<spaceDim;j++){
    if(max[j]>min[j]) boxDim=j+1;
  }
}