polyhedron.h

another source code for reduced basis

/*
A fast and simple stretch-minimizing mesh parameterization C++ code
Copyright:(c) Shin Yoshizawa, 2004
E-mail: [e-mail]
URL: http://www.mpi-sb.mpg.de/~shin
Affiliation: Max-Planck-Institut fuer Informatik: Computer Graphics Group 
 Stuhlsatzenhausweg 85, 66123 Saarbruecken, Germany
 Phone [phone] Fax [phone] 

 All right is reserved by Shin Yoshizawa.
This C++ sources are allowed for only primary user of 
research and educational purposes. Don't use secondary: copy, distribution, 
diversion, business purpose, and etc.. 
 */
#define PI 3.1415926535897932385
#define next(p) p->next
#define back(p) p->back
// parameterization class
class Polyhedron{
 public:
  int numberV;
  int numberF;
  int boundarytype;
  /* 
     0: unit square
     1: unit circle
     2: natural boundary of the intrinsic parameterization
  */
  int paramtype;
  /* 
     1: U1 parameterization: single optimization step
     2: Uopt parameterization: optimal parameterization
  */
  int weighttype;
  /*
    0: Shape Preserving: Floater CAGD 1997.
    1: w = 1: Tutte 1963. 
    2: Harmonic Map: Eck et al. SIGGRAPH 1995.
    3: Intrinsic Map: Desbrun et al. EUROGRAPHICS 2002.
       ((1-intrinsiclambda)*Chi + intrinsiclambda*Authalic)
       If intrinsiclambda=1 then it is equivalent to Eck's map
    4: Mean value: Floater CAGD 2003.

   */
  double intrinsiclambda;
  
  int smooth;
  /* 
     For the optimal low stretch parameterization,
     1: fast and aggressive
     2: smooth: if you will get a problem (unstability) using 1 then it is better to use 2.
  */  
  Point3d **point;
  int pickID;/* first boundary point ID: it will be conner of the unit square. 
		if pickID=-1 then pickID will be the first boundary vertex ID.
	      */

  int boundarysigma;/* When we test current stretch distortion,
		       0 without boundary triangles,
		       1 with boundary triangles
		    */
  int level;  

  /* Parameter Mesh Coordinates */
  double *pU;
  double *pV;
  /*******************************/

  int **Face;
  double *areaMap3D;
  double sumarea3D;
  double constsumarea3D;
  
  double *sigma;
  double *E;
  double *G;
  
  IDSet *IDtool;
  PointTool *PT;
  /*************************/
  
  // 1-ring neighbor structure
  PolarList **PHead;
  PolarList **PTail;
  IDList **IHead;
  IDList **ITail;
  IDList **VHead;
  IDList **VTail;
  IDList **FHead;
  IDList **FTail;
  int *boundary;
  int *neighborI;
  int *neighborF;
  
  int numboundary;
  Point3d **bc;
  Point3d **dbc;  

  double totalarea3D;
  double totalarea2D;
  double gammaP;/* transfer parameter: eta in the paper */

  int iteNum; /* Maximum iteration number in PCBCG */
  double PCBCGerror;/* Error threshold in PCBCG */
  
  
  
  Polyhedron(){
    // Default setting
    pickID=-1;
    paramtype=2;
    boundarytype=0;
    weighttype=0;
    iteNum=2000;
    gammaP=1.0;
    PCBCGerror=pow(0.1,6.0);
    smooth=1;
    intrinsiclambda=0.5;
    boundarysigma=1;
  }
  virtual ~Polyhedron(){
    if(numberV!=0&&point!=NULL){
      memorydelete();
    }
  }
  // parameterization
  void Polyhedron::param(){
    //boundary mapping 
    
    BoundaryMap();
    
    if(boundarytype==2){
      setNaturalB(iteNum,PCBCGerror); 
      //make local parameterization via polar map 
      setPolarMap(); 
    }else{
      //make local parameterization via polar map 
      setPolarMap(); 
    }
    
    if(paramtype==1){
      ParametrizationSingle(iteNum,PCBCGerror);
      

    }else{
      if(smooth==1){ 
	ParametrizationOptimal(iteNum,PCBCGerror);
      }else{ 
	ParametrizationSmoothOptimal(iteNum,PCBCGerror);
      }
    }
     
  }

  
  void Polyhedron::readmesh(char *filename){
    
    FILE *in=NULL;
    in = fopen(filename,"r");
    int dV=0;
    int dF=0;
    int i,j;
    int di=0;
    int dj=0;
    int dk=0;
    double dx=0.0;
    double dy=0.0;
    double dz=0.0;
    fscanf(in,"%d",&dV);
    fscanf(in,"%d",&dF);
    
    memoryallocate(dV,dF);
    for(i=0;i<numberV;i++){
      fscanf(in,"%lf %lf %lf",&dx,&dy,&dz);
      setPoint(i,dx,dy,dz);
    }
    int val=3;
    for(i=0;i<numberF;i++){
      fscanf(in,"%d %d %d %d",&val,&di,&dj,&dk);
      setFace(i,di,dj,dk);
      IDtool->AppendVFSort(i,FHead[Face[i][0]],FTail[Face[i][0]]);
      IDtool->AppendVFSort(i,FHead[Face[i][1]],FTail[Face[i][1]]);
      IDtool->AppendVFSort(i,FHead[Face[i][2]],FTail[Face[i][2]]);
    }
    
    
    
    fclose(in);
    
    /* feature analysis */
    
    SetBoundaryLines();
    setAreaMap3D();
  }
  void Polyhedron::writemesh(char *filename){
    int i=0;
    FILE *out = fopen(filename,"w");
    fprintf(out,"%d\n",numberV);
    fprintf(out,"%d\n",numberF);
    for(i=0;i<numberV;i++){
      fprintf(out,"%lf %lf %lf\n",pU[i],pV[i],0.0);
    }
    for(i=0;i<numberF;i++)
      fprintf(out,"3 %d %d %d\n",Face[i][0],Face[i][1],Face[i][2]);
    fclose(out);
    
  }
  

  Polyhedron(const Polyhedron& rhs);
  const Polyhedron &operator=(const Polyhedron& rhs);




private:



  void Polyhedron::BoundaryMap(){
    int i,j,k;
  
    pU = new double[numberV];
    pV = new double[numberV];
    double maxvalU,minvalU,valrangeU;
    double maxvalV,minvalV,valrangeV;
    
    IDList *now=NULL;
    IDList *BedgeH = new IDList();
    IDList *BedgeT = new IDList();
    BedgeH->next = BedgeT;
    BedgeT->back = BedgeH;
    int cnt=0;
    
    for(i=0;i<numberF;i++){
      if(boundary[Face[i][0]]==1&&
	 boundary[Face[i][1]]==1&&
	 boundary[Face[i][2]]==1){
	
	  IDtool->AppendVF(Face[i][0],BedgeT);
	  IDtool->AppendVF(Face[i][1],BedgeT);
	  IDtool->AppendVF(Face[i][1],BedgeT);
	  IDtool->AppendVF(Face[i][2],BedgeT);
	  IDtool->AppendVF(Face[i][2],BedgeT);
	  IDtool->AppendVF(Face[i][0],BedgeT);
	
	
	
      }else{
	if((boundary[Face[i][0]]==1&&boundary[Face[i][1]]==1)||
	   (boundary[Face[i][1]]==1&&boundary[Face[i][2]]==1)||
	   (boundary[Face[i][2]]==1&&boundary[Face[i][0]]==1)
	   ){
	  if(boundary[Face[i][0]]==1&&boundary[Face[i][1]]==1){
	    IDtool->AppendVF(Face[i][0],BedgeT);
	    IDtool->AppendVF(Face[i][1],BedgeT);
	    
	  }else if(boundary[Face[i][1]]==1&&boundary[Face[i][2]]==1){
	    IDtool->AppendVF(Face[i][1],BedgeT);
	    IDtool->AppendVF(Face[i][2],BedgeT);
	    
	  }else if(boundary[Face[i][2]]==1&&boundary[Face[i][0]]==1){
	    IDtool->AppendVF(Face[i][2],BedgeT);
	    IDtool->AppendVF(Face[i][0],BedgeT);
	    
	  }
	}
      }
    }
    
    int jumpcheck=0;
    now = BedgeH;
    while(next(now)!=BedgeT){
      now = next(now);jumpcheck=0;
      if(boundary[now->ID]== 1&&boundary[next(now)->ID]==1){
	int checkbedge = 0;
	IDList *now2 = VHead[now->ID];
	while(next(now2)!=VTail[now->ID]){
	  now2 = next(now2);
	  if(now2->ID==next(now)->ID||
	     next(now2)->ID==next(now)->ID){
	    checkbedge++;
	  }
	  now2 = next(now2);
	  
	}
	if(checkbedge>=2){
	  
	  IDList *dummyn = next(now)->next;
	  IDList *dummyb = now->back;
	  dummyb->next = dummyn;
	  dummyn->back = dummyb;
	  
	  IDList *ddn = next(now);
	  IDList *ddnn = now;
	  delete ddn;
	  delete ddnn;
	  
	  now = dummyb;
	  jumpcheck=1;
	}
      }
      if(jumpcheck!=1)
      now = next(now);
    }
    

    
    
    IDList *BpointH = new IDList();
    IDList *BpointT = new IDList();
    
    BpointH->next = BpointT;
    BpointT->back = BpointH;
    double tlength=0.0;
    now = BedgeH;
    while(next(now)!=BedgeT){
      now = next(now);
      tlength += PT->Distance(point[now->ID],point[next(now)->ID]);    
      
      now = next(now);
      
    }
    int startID=0;
    int *checkbin = new int[numberV];
    for(i=0;i<numberV;i++)checkbin[i]=-1;
    if(pickID==-1){
      for(i=0;i<numberV;i++)if(boundary[i]==1){
	pickID=i;break;
      }
    }
    
    //printf("start id = %d\n",pickID);
    startID=pickID;
    checkbin[pickID]=0;
    
  IDtool->AppendVF(pickID,BpointT);
  k=1;
  
  
  while(1){
    now = BedgeH;
    
    while(next(now)!=BedgeT){
      now = next(now);
      if(startID==now->ID||startID==next(now)->ID){
	if(startID==now->ID&&checkbin[next(now)->ID]==-1){
	  k++;
	  checkbin[next(now)->ID]=0;
	  IDtool->AppendVF(next(now)->ID,BpointT);
	  startID=next(now)->ID;
	  break;
	}else if(startID==next(now)->ID&&checkbin[now->ID]==-1){
	  k++;
	  checkbin[now->ID]=0;
	  IDtool->AppendVF(now->ID,BpointT);
	  startID=now->ID;
	  break;
	}
	
      }
      now = next(now);
    }
    
    if(k>=numboundary)break;
  }
  
  IDtool->CleanNeighbor(BedgeH,BedgeT);
  tlength *= 0.25;
  
  double clen=0.0;
  int state=0;
  i=0;
  
  
  /* Unit Circle */
  if(boundarytype==1){
    now = BpointH;
    while(next(now)!=BpointT){
      now = next(now);
      pU[now->ID] = 0.5*cos(2.0*PI*(((double)(i))/((double)(numboundary))));
      pV[now->ID] = 0.5*sin(2.0*PI*(((double)(i))/((double)(numboundary))));
      i++;
    }
  }else{
    
    /* Unit Square */
    
    now = BpointH;
    
    clen=0.0;
    
    
    while(next(next(now))!=BpointT){
      now = next(now);
      
      if(clen<1.0){
	pU[now->ID] = clen; 
	pV[now->ID]=0.0;
	
      }else if(clen>=1.0&&clen<2.0){
	if(state==0){
	  clen = 1.0;
	  pU[now->ID] = 1.0;
	  pV[now->ID]=0.0;
	  
	  state=1;
	}else{
	  pU[now->ID] = 1.0;
	  pV[now->ID]= clen-1.0;
	  
	}
	
      }else if(clen>=2.0&&clen<3.0){
	if(state==1){
	  clen = 2.0;
	  pU[now->ID] = 1.0;
	  pV[now->ID] = 1.0;
	  state=2;
	}else{
	  pU[now->ID] = 3.0-clen;
	  pV[now->ID]= 1.0;
	}
	
	
      }else if(clen>=3.0){
	if(state==2){
	  clen = 3.0;
	  pU[now->ID] = 0.0;
	  pV[now->ID] = 1.0;
	  state=3;
	}else{
	  pU[now->ID] = 0.0;
	  pV[now->ID] = 4.0-clen;
	}
      }
      
      clen += PT->Distance(point[now->ID],point[next(now)->ID])/tlength;
    }
    now = next(now);
    clen += PT->Distance(point[now->ID],point[next(BpointH)->ID])/tlength;
    pU[now->ID] = 0.0;
    pV[now->ID] = 4.0-clen;
  }
  
  for(i=0;i<numberV;i++){
    if(boundary[i]==1){
            
    }else{
      pU[i] *= 0.5;
      pV[i] *= 0.5;    
    }
  }
  
  IDtool->CleanNeighbor(BpointH,BpointT);
  delete [] checkbin;
  
  
  
  
  
  }
  
  void Polyhedron::memorydelete(){
    int i;
    if(E!=NULL)delete [] E;
    if(G!=NULL)delete [] G;
    for(i=0;i<10;i++){
      delete bc[i];
      delete dbc[i];
    }
    delete [] bc;
    bc = NULL;
    delete [] dbc;
    dbc = NULL;
    
    if(pU!=NULL){
      delete [] pU;
      pU = NULL;
      
    }
    if(pV!=NULL){
      delete [] pV;
      pV = NULL;
      
    }
    if(areaMap3D!=NULL){
      delete [] areaMap3D;
      
    }
    if(IDtool!=NULL){
      if(FHead!=NULL){
	IDtool->CleanNeighborL(FHead,FTail,numberV);
	
      }
      if(IHead!=NULL&&ITail!=NULL){
	IDtool->CleanNeighborL(IHead,ITail,numberV);
	IHead=NULL;
	ITail=NULL;
      }
      if(PHead!=NULL&&PTail!=NULL){
	for(i=0;i<numberV;i++){
	  
	  if(boundary[i]!=1&&PHead[i]!=NULL&&PTail[i]!=NULL&&PHead[i]->next!=PTail[i]){
	    IDtool->CleanNeighborPolar(PHead[i],PTail[i]);
	  }
	}
	delete [] PHead;
	delete [] PTail;
	
      }
      
      if(VHead!=NULL&&VTail!=NULL){
	IDtool->CleanNeighborL(VHead,VTail,numberV);
	VHead=NULL;
	VTail=NULL;
      }
      delete IDtool;
      IDtool=NULL; 
    }
    
  if(point!=NULL){
    if(numberV!=0){
      for(i=0;i<numberV;i++){
	delete point[i];
      }
    }
    delete [] sigma;
    delete [] point;
    point = NULL;
    
  }
  if(Face!=NULL){
    if(numberF!=0){
      for(i=0;i<numberF;i++)delete [] Face[i];
    }
    delete [] Face;
    Face=NULL;
  }
  if(neighborI!=NULL){
    delete [] neighborI;
  }
 if(neighborF!=NULL){
    delete [] neighborF;
  }
 if(boundary!=NULL){
    delete [] boundary;
  }
 delete PT;
  }  
  
  void Polyhedron::memoryallocate(int dV,int dF){
            
    numberV=dV;
    numberF=dF;
    E = new double[numberF];
    G = new double[numberF];
    sigma = new double[numberV];
    point = new Point3d* [numberV];
    Face = new int* [numberF];
    IDtool = new IDSet();
    IHead = new IDList* [numberV];
    ITail = new IDList* [numberV];
    PHead = new PolarList* [numberV];
    PTail = new PolarList* [numberV];
  
    VHead = new IDList* [numberV];
    VTail = new IDList* [numberV];
    FHead= new IDList* [numberV];
    FTail= new IDList* [numberV];
    neighborF= new int[numberV];
    boundary = new int[numberV];
    neighborI = new int[numberV]; 
    areaMap3D = new double[numberF];
    
    numboundary=0;
    bc = new Point3d* [10];
    dbc = new Point3d* [10];
    int i;
    for(i=0;i<10;i++){
      bc[i] = new Point3d(0.0,0.0,0.0);
      dbc[i] = new Point3d(0.0,0.0,0.0);
    }
  }
  void Polyhedron::setPoint(int i,double dx,double dy,double dz){
  point[i] = new Point3d(dx,dy,dz);
  IHead[i] = new IDList();