DELAUNAY.C

Go to the documentation of this file.
00001 /* --
00002 OpenFX version 2.0 - Modelling, Animation and Rendering Package
00003 Copyright (C) 2000 - 2007 OpenFX Development Team
00004 -- */
00005 
00006 /*  DELAUNAY.C  */
00007 
00008 #include <stdlib.h>
00009 #include <stdio.h>
00010 #include <math.h>
00011 #include <windows.h>
00012 
00013 #include "struct.h"
00014 #include "dstruct.h"
00015 
00016 #define TRITOP           0
00017 #define TRIFRONT         1
00018 #define TRIRIGHT         2
00019 #define MAXUNIT    4194304L     //1073741824L
00020 
00021 #define DLG_CUT                     600
00022 #define DLG_CUT_SPLIT_Y             601
00023 #define DLG_CUT_SPLIT_N             602
00024 #define DLG_CUT_BOTH                603
00025 #define DLG_CUT_WORKPIECE           604
00026 
00027 
00028 #include "delaunay.h"
00029 
00030 #define DESELECTED 0
00031 #define SELECTED   1
00032 #define INEDITOR   3
00033 
00034 static HWND      hParent;
00035 static HINSTANCE hThisInstance;
00036 
00037 
00038 #if __WATCOMC__
00039 int APIENTRY LibMain(HANDLE hDLL, DWORD dwReason, LPVOID lpReserved){
00040 #else
00041 BOOL WINAPI DllMain(HANDLE hDLL, DWORD dwReason, LPVOID lpReserved){
00042 #endif
00043   switch (dwReason) {
00044     case DLL_PROCESS_ATTACH: {
00045       hThisInstance=hDLL;
00046       if(hDLL == NULL)MessageBox( GetFocus()," NULL instance",NULL, MB_OK);
00047       break;
00048     }
00049     case DLL_PROCESS_DETACH:
00050       break;
00051   }
00052   return TRUE;
00053 }
00054 
00055 #if __SC__
00056 #pragma startaddress(DllMain)
00057 #endif
00058 
00059 static BOOL CALLBACK BooleanDlgProc(HWND hwnd, UINT msg,
00060                              WPARAM wparam, LPARAM lparam){
00061  BOOL err;
00062  int i;
00063  static struct bp {BOOL sp,al;} *dp;
00064  switch( msg ) {
00065    case WM_INITDIALOG:
00066      dp=(struct bp *)lparam;
00067      CentreDlgOnS(hwnd);
00068      return (TRUE);
00069    case WM_COMMAND:
00070       switch(LOWORD(wparam)){
00071         case IDCANCEL:
00072           EndDialog(hwnd,FALSE);
00073           return(TRUE);
00074         case IDOK:
00075           EndDialog(hwnd,TRUE);
00076           return(TRUE);
00077         default:
00078           break;
00079       }
00080       break;
00081     default: break;
00082  }
00083  return(FALSE);
00084 }
00085 
00086 BOOL _Xmodeler
00087  (HWND parent_window,HWND info_window,X__STRUCTURE *lpevi){
00088  HCURSOR hSave;
00089  int type,split;
00090  struct bp {BOOL sp,al;} d;
00091  lpEVI=lpevi;
00092  hParent=parent_window;
00093  d.sp=TRUE; d.al=TRUE;
00094  if(DialogBoxParam(hThisInstance,MAKEINTRESOURCE(DLG_CUT),hParent,
00095              (DLGPROC)BooleanDlgProc,(LPARAM)&d) == FALSE)return FALSE;
00096  if(d.al)type=0;  else type=1;
00097  if(d.sp)split=1; else split=0;
00098  hSave=SetCursor(LoadCursor(NULL,IDC_WAIT));
00099  if     (ActiveView == TRITOP)IntersectObject(2);
00100  else if(ActiveView == TRIFRONT)IntersectObject(1);
00101  else if(ActiveView == TRIRIGHT)IntersectObject(0);
00102  SetCursor(hSave);
00103  return TRUE;
00104 }
00105 
00106 #define TOL  200000.0  // very important DO not TAMPER
00107 #define TOL1    0.0005 //0.001 // 0.01 changed 4/6/95 for problem with missed point
00108 #define TOL2    0.9995 //0.999 // 0.99 in phase one too near an edge and disregarded
00109 #define TOL3   -1.00000
00110 #define TOL8    1.00000
00111 #define TOL5   -0.001
00112 #define TOL6    1.001
00113 #define TOL9  100.0
00114 #define TOLA  100.0
00115 #define TOLB  10.0
00116 #define TOLC  100.0
00117 
00118 static short EdgePierceFace(vector vi,
00119                             vector p0, vector p1, vector p2,
00120                             short *code, short *id){
00121 /* code  = 0 point vi is in plane and inside face                */
00122 /*       = 1   "      "      "        on edge id id's the edge   */
00123 /*       = 2   "      "      "        at vertex id id's the edge */
00124  double ve1,ve2,ve3,vp1,vp2,vp3,q1,q2,q3,det,a,b,ab,
00125         det1,det2,det3,detmax,dm;
00126  short k;
00127   q1=vi[0]-p0[0];  q2=vi[1]-p0[1];  q3=vi[2]-p0[2];
00128  ve1=p2[0]-p0[0]; ve2=p2[1]-p0[1]; ve3=p2[2]-p0[2];
00129  vp1=p1[0]-p0[0]; vp2=p1[1]-p0[1]; vp3=p1[2]-p0[2];
00130  det1=ve1*vp2-vp1*ve2;    /* new way to ensure optimum numerics */
00131  det2=ve2*vp3-vp2*ve3;
00132  det3=ve1*vp3-vp1*ve3;
00133  k=0; detmax=TOL;
00134  if((dm=fabs(det1)) > detmax){k=1; detmax=dm;}
00135  if((dm=fabs(det2)) > detmax){k=2; detmax=dm;}
00136  if((dm=fabs(det3)) > detmax){k=3; detmax=dm;}
00137  if(k == 0)return 0;
00138  else if(k == 1){
00139    a=( vp2*q1-vp1*q2)/det1;
00140    b=(-ve2*q1+ve1*q2)/det1;
00141  }
00142  else if(k == 2){
00143    a=( vp3*q2-vp2*q3)/det2;
00144    b=(-ve3*q2+ve2*q3)/det2;
00145  }
00146  else if(k == 3){
00147    a=( vp3*q1-vp1*q3)/det3;
00148    b=(-ve3*q1+ve1*q3)/det3;
00149  }
00150  else return 0;
00151 /* a=1 at p2      b=1 at p1  */
00152  ab=a+b;
00153  if(a < TOL5 || a > TOL6 || b < TOL5 || b > TOL6 || ab > TOL6)return 0;
00154  if(a >= TOL1 && a <= TOL2 && b >= TOL1 && b <= TOL2 && ab <= TOL2){
00155    *code=0;  /* inside face */
00156    return 1;
00157  }
00158  else if(a < TOL1){ /* along edge 0-1 */
00159    if(b < TOL1){ /* p0 */
00160      *code=2; *id=0;
00161      return 1;
00162    }
00163    else if(b > TOL2){ /* p1 */
00164      *code=2; *id=1;
00165      return 1;
00166    }
00167    else{  /* on edge 0-1 */
00168      *code=1; *id=0;
00169      return 1;
00170    }
00171  }
00172  else if(b < TOL1){ /* along edge 0-2 */
00173    if(a < TOL1){ /* p0 */
00174      *code=2; *id=0;
00175      return 1;
00176    }
00177    else if(a > TOL2){ /* p2 */
00178      *code=2; *id=2;
00179      return 1;
00180    }
00181    else{ /* on edge 0-2 */
00182      *code=1; *id=2;
00183      return 1;
00184    }
00185  }
00186  else if(ab > TOL2){/* on edge 1-2 */
00187    *code=1; *id=1;
00188    return 1;
00189  }
00190  return 0;
00191 }
00192 
00193 static int O_normalize(vector v){
00194  double n,nn;
00195  n= (double)((v[0]*v[0]) + (v[1]*v[1]) + (v[2]*v[2]));
00196  if(n < 1.e-10)return FAIL;
00197  nn=sqrt(n);
00198  n = 1.0 / nn;
00199  VECSCALE(n,v,v)
00200  return OK;
00201 }
00202 
00203 static short Del_Swap(vector p, long v1l, long v2l, long v3l){
00204  /* test to see if edge common to triangle pair  v1,v2,v3  and v1,v2,p */
00205  /* should be swapped so as to keep the triangles as near equilateral  */
00206  /* as possible.                                                       */
00207  vertex *v1,*v2,*v3;
00208  double cosa,cosb,sina,sinb,sinab;
00209  vector vp3,v13,v23,va,vb,v1p,v2p;
00210  v1=(MainVp+v1l); v2=(MainVp+v2l); v3=(MainVp+v3l);
00211  VECSUB((double)v1->xyz,(double)v3->xyz,v13) O_normalize(v13);
00212  VECSUB((double)v2->xyz,(double)v3->xyz,v23) O_normalize(v23);
00213  VECSUB((double)v1->xyz,p,v1p) O_normalize(v1p);
00214  VECSUB((double)v2->xyz,p,v2p) O_normalize(v2p);
00215  cosa=DOT(v13,v23);
00216  cosb=DOT(v1p,v2p);
00217  if(cosa >= 0.0 && cosb >= 0.0)return 0;     /* a < 90 and b < 90 */
00218  if(cosa <  0.0 && cosb <  0.0)return 1;
00219  CROSS(v13,v23,va) sina=va[0]*va[0]+va[1]*va[1]+va[2]*va[2];
00220  if(sina < 1.e-10)return 0; sina = sqrt(sina);
00221  CROSS(v2p,v1p,vb) sinb=vb[0]*vb[0]+vb[1]*vb[1]+vb[2]*vb[2];
00222  if(sinb < 1.e-10)return 0; sinb = sqrt(sinb);
00223  sinab=sina*cosb+sinb*cosa;
00224  if(sinab < 0.0)return 1;
00225  return 0;
00226 }
00227 
00228 static long BoolAddVertex(vector v, short status){
00229  vertex *vv;
00230  if(!UpdateVertexHeap(Nvert+1)){
00231    MessageBeep(MB_OK);
00232    return Nvert-1;
00233  }
00234  CreateVertex();
00235  vv=(MainVp+Nvert-1);
00236  vv->xyz[0]=(long)v[0];
00237  vv->xyz[1]=(long)v[1];
00238  vv->xyz[2]=(long)v[2];
00239  vv->id=0; vv->status=status;
00240  if(status == DESELECTED){NvertSelect--; NvertDeselect++;}
00241  return Nvert-1;
00242 }
00243 
00244 static FACET *add_small_face(long a, long b, long c,
00245                              long f, long i, FACET *w,
00246                              long *nface){
00247   long nf;
00248   nf = *nface;
00249   if((w=(FACET *)X__Realloc(w,(nf+1)*(long)sizeof(FACET)))
00250        == NULL)return NULL;
00251   if(!UpdateFaceHeap(Nface+1))return NULL;
00252   CreateFace(a,b,c);
00253   CopyFaceProp((MainFp+f),(MainFp+Nface-1));
00254   w[nf].f=Nface-1;
00255   (*nface)++;
00256   return w;
00257 }
00258 
00259 static long Del_Edgeid(long v1, long v2, FACET *w, long k,
00260                        long *fa1, long *fa2, long *e1, long *e2){
00261  face *f;
00262  f=(MainFp+w[k].f);
00263  if((f->V[0] == v1 && f->V[1] == v2) ||
00264     (f->V[0] == v2 && f->V[1] == v1)){
00265    *fa1=w[k].fa[1]; *fa2=w[k].fa[2];
00266    *e1=w[k].ea[1]; *e2=w[k].ea[2];
00267    return f->V[2];
00268  }
00269  else if((f->V[1] == v1 && f->V[2] == v2) ||
00270          (f->V[1] == v2 && f->V[2] == v1)){
00271    *fa1=w[k].fa[2]; *fa2=w[k].fa[0];
00272    *e1=w[k].ea[2]; *e2=w[k].ea[0];
00273    return f->V[0];
00274  }
00275  else if((f->V[2] == v1 && f->V[0] == v2) ||
00276          (f->V[2] == v2 && f->V[0] == v1)){
00277    *fa1=w[k].fa[0]; *fa2=w[k].fa[1];
00278    *e1=w[k].ea[0]; *e2=w[k].ea[1];
00279    return f->V[1];
00280  }
00281  MessageBox(NULL,"Algorithm confused",NULL,MB_OK);
00282  return -1;
00283 }
00284 
00285 static FACET *AddDelaunayTriangulation(vector p2,long *nf, FACET *w){
00286  short code,id;
00287  long j,k,l,m,t,n,nn,*Dstack=NULL,Dss=0,fa1,fa2;
00288  FACET *bf;
00289  vector pa,pb,pc;
00290  long vp,w0,w1,w2,vo,vt;
00291  long f1;
00292  long e0,e1,e2,ep;
00293  n= *nf;
00294  for(j=0;j<n;j++){ /* find which face this intersects    */
00295    bf=(w+j);
00296    f1=bf->f;
00297    w0=(MainFp+f1)->V[0]; w1=(MainFp+f1)->V[1]; w2=(MainFp+f1)->V[2];
00298    POINT2VECTOR((MainVp+w0)->xyz,pa)
00299    POINT2VECTOR((MainVp+w1)->xyz,pb)
00300    POINT2VECTOR((MainVp+w2)->xyz,pc)
00301    if(EdgePierceFace(p2,pa,pb,pc,&code,&id) /* && code == 0 */){
00302      vp=BoolAddVertex(p2,DESELECTED);
00303      if(!UpdateEdgeHeap(Nedge+3))return NULL;
00304      CreateEdge(vp,w0);  e0=Nedge-1;
00305      CreateEdge(vp,w1);  e1=Nedge-1;
00306      CreateEdge(vp,w2);  e2=Nedge-1;
00307      if((w=add_small_face(vp,w0,w1,f1,j,w,nf))==NULL)return NULL;
00308      if((w=add_small_face(vp,w2,w0,f1,j,w,nf))==NULL)return NULL;
00309      (MainFp+f1)->V[0]=vp;       /* modify this face       */
00310      nn = *nf;
00311      w[nn-2].fa[0]=nn-1;         /* adjacency in new faces */
00312      w[nn-2].fa[1]=w[j].fa[0];
00313      w[nn-2].fa[2]=j;
00314      w[nn-1].fa[0]=j;
00315      w[nn-1].fa[1]=w[j].fa[2];
00316      w[nn-1].fa[2]=nn-2;
00317      if((k=w[j].fa[0]) >= 0){   /* change adjacency in adjacent */
00318        for(l=0;l<3;l++)if(w[k].fa[l] == j)w[k].fa[l]=nn-2;
00319      }
00320      if((k=w[j].fa[2]) >= 0){
00321        for(l=0;l<3;l++)if(w[k].fa[l] == j)w[k].fa[l]=nn-1;
00322      }
00323      w[j].fa[0]=nn-2;           /*  change adjacency in this face */
00324      w[j].fa[2]=nn-1;
00325      w[nn-2].ea[0]=e0;
00326      w[nn-2].ea[1]=w[j].ea[0];
00327      w[nn-2].ea[2]=e1;
00328      w[nn-1].ea[0]=e2;
00329      w[nn-1].ea[1]=w[j].ea[2];
00330      w[nn-1].ea[2]=e0;
00331      w[j].ea[0]=e1;
00332      w[j].ea[2]=e2;
00333      if(Dstack == NULL){
00334        if((Dstack=(long *)X__Malloc(sizeof(long)*(3)))
00335                  == NULL)return NULL;
00336      }
00337      else {
00338        if((Dstack=(long *)X__Realloc(Dstack,sizeof(long)*(Dss+3)))
00339                  == NULL)return NULL;
00340      }
00341      if(w[j].fa[1] >= 0)   { Dstack[Dss] = j;    Dss++;}
00342      if(w[nn-2].fa[1] >= 0){ Dstack[Dss] = nn-2; Dss++;}
00343      if(w[nn-1].fa[1] >= 0){ Dstack[Dss] = nn-1; Dss++;}
00344      while(Dss > 0){
00345        Dss--;
00346        j=Dstack[Dss]; k=w[j].fa[1];
00347        if(Dss > 0){
00348          if((Dstack=(long *)X__Realloc(Dstack,sizeof(long)*(Dss))) == NULL){
00349            return NULL;  /* pop off top of stack */
00350          }
00351        }
00352        if((vo=Del_Edgeid((MainFp+(w[j].f))->V[1],
00353                          (MainFp+(w[j].f))->V[2],w,k,
00354                          &fa1,&fa2,&e1,&e2)) < 0)return NULL;
00355        if(Del_Swap(p2,(MainFp+(w[j].f))->V[1],(MainFp+(w[j].f))->V[2],vo)){
00356          e0=w[j].ea[1];
00357          if((m=w[j].fa[2]) >= 0){   /* change adjacency in adjacent */
00358            for(l=0;l<3;l++)if(w[m].fa[l] == j)w[m].fa[l]=k;
00359          }
00360          if(fa1 >= 0){
00361            for(l=0;l<3;l++)if(w[fa1].fa[l] == k)w[fa1].fa[l]=j;
00362          }
00363          (MainEp+e0)->V[0]=vp;
00364          (MainEp+e0)->V[1]=vo; /* reverse edge */
00365          vt=(MainFp+(w[j].f))->V[2];
00366          (MainFp+(w[j].f))->V[2]=vo;   /* update face j */
00367          w[j].fa[1]=fa1; t=w[j].fa[2]; w[j].fa[2]=k;
00368          w[j].ea[1]=e1; ep=w[j].ea[2]; w[j].ea[2]=e0;
00369          (MainFp+(w[k].f))->V[0]=vp;
00370          (MainFp+(w[k].f))->V[1]=vo;
00371          (MainFp+(w[k].f))->V[2]=vt; /* face k */
00372          w[k].fa[0]=j; w[k].fa[1]=fa2; w[k].fa[2]=t;
00373          w[k].ea[0]=e0; w[k].ea[1]=e2; w[k].ea[2]=ep;
00374          if((Dstack=(long *)X__Realloc(Dstack,sizeof(long)*(Dss+2))) == NULL){
00375            return NULL;  /* push both onto stackstack */
00376          }
00377          if(w[j].fa[1] >= 0){Dstack[Dss] = j; Dss++;}
00378          if(w[k].fa[1] >= 0){Dstack[Dss] = k; Dss++;}
00379        }
00380      }
00381      if(Dstack != NULL)X__Free(Dstack);
00382      return w; /* intersection found */
00383    }
00384  }
00385  return w;
00386 }
00387 
00388 static void IntersectObject(int plane){
00389  /* FACET is the face being worked on           */
00390  vector *Points=NULL;
00391  FACET *w=NULL;
00392  vertex *vpp;
00393  long w0,w1,w2;
00394  long i,j,k,l,m;
00395  vector p0,p1,p2,du,n;
00396  long n_work_face=0,nPoints=0;
00397  long xc, d_max_x = -MAXUNIT, d_min_x = MAXUNIT;
00398  long yc, d_max_y = -MAXUNIT, d_min_y = MAXUNIT;
00399  long zc, d_max_z = -MAXUNIT, d_min_z = MAXUNIT;
00400  long radius;
00401  if(NvertSelect   == 0)return;
00402  if(NvertDeselect != 0){
00403    vpp=MainVp; for(i=0;i<Nvert;i++){
00404      if(vpp->status == DESELECTED){
00405        vpp->status=INEDITOR;
00406      }
00407      vpp++;
00408    }
00409  }
00410  //
00411  // Make the bounding first triangle
00412  //
00413  l=plane;
00414  vpp=MainVp; for(i=0;i<Nvert;i++){
00415    if(vpp->status == SELECTED){
00416      if(l == 2){j=0; k=1; m=2;}
00417      if(l == 0){j=1; k=2; m=0;}
00418      if(l == 1){j=2; k=0; m=1;}
00419      d_min_x=min(d_min_x,vpp->xyz[j]);
00420      d_max_x=max(d_max_x,vpp->xyz[j]);
00421      d_min_y=min(d_min_y,vpp->xyz[k]);
00422      d_max_y=max(d_max_y,vpp->xyz[k]);
00423      d_min_z=min(d_min_z,vpp->xyz[m]);
00424      d_max_z=max(d_max_z,vpp->xyz[m]);
00425    }
00426    vpp->id=0; vpp++;
00427  }
00428  xc = (d_min_x/2 + d_max_x/2);
00429  yc = (d_min_y/2 + d_max_y/2);
00430  zc = (d_min_z/2 + d_max_z/2);
00431  radius = max((d_max_x - d_min_x),(d_max_y - d_min_y));
00432  n[0]=n[1]=n[2]=0.0;
00433  if(l == 2){i=0; j=1; k=2; n[2]=1.0;}
00434  if(l == 0){i=1; j=2; k=0; n[0]=1.0;}
00435  if(l == 1){i=2; j=0; k=1; n[1]=1.0;}
00436  p0[i]=(double)xc;
00437  p0[j]=(double)yc+2.0*(double)radius;
00438  p0[k]=zc;
00439  p1[i]=(double)xc-1.732*(double)radius;
00440  p1[j]=(double)yc-(double)radius;
00441  p1[k]=zc;
00442  p2[i]=(double)xc+1.732*(double)radius;
00443  p2[j]=(double)yc-(double)radius;
00444  p2[k]=zc;
00445 
00446  BoolAddVertex(p0,DESELECTED); w0=Nvert-1;
00447  BoolAddVertex(p1,DESELECTED); w1=Nvert-1;
00448  BoolAddVertex(p2,DESELECTED); w2=Nvert-1;
00449  if(!UpdateEdgeHeap(Nedge+3))goto QUIT;
00450  CreateEdge(w0,w1);
00451  CreateEdge(w1,w2);
00452  CreateEdge(w2,w0);
00453  if(!UpdateFaceHeap(Nface+1))goto QUIT;
00454  CreateFace(w0,w1,w2);
00455  if((w=(FACET *)X__Malloc((long)sizeof(FACET)))
00456        == NULL)goto BERR;
00457  n_work_face=0;
00458  w[n_work_face].f=Nface-1;
00459  for(i=0;i<3;i++){
00460    w[n_work_face].fa[i] = -1;
00461    w[n_work_face].ea[i] =  i;
00462  }
00463  n_work_face++;
00464  if((Points=(vector *)X__Malloc(NvertSelect*(long)sizeof(vector)))
00465       == NULL)goto BERR;
00466  nPoints=0;
00467  for(vpp=MainVp,i=0;i<Nvert;i++){
00468    if(vpp->status == SELECTED){
00469      POINT2VECTOR(vpp->xyz,Points[nPoints])
00470      nPoints++;
00471    }
00472    vpp++;
00473  }
00474  for(i=0;i<nPoints;i++){
00475    double mu;
00476    VECSUB(Points[i],p0,du)
00477    mu=DOT(du,n);
00478    VECSCALE(mu,n,du)
00479    VECSUM(Points[i],du,p2)
00480    if((w=AddDelaunayTriangulation(p2,&n_work_face,w)) == NULL)goto BERR;
00481  }
00482  (MainVp+w0)->status=SELECTED;
00483  (MainVp+w1)->status=SELECTED;
00484  (MainVp+w2)->status=SELECTED;
00485  NvertSelect   += 3;
00486  NvertDeselect -= 3;
00487  EraseVertex(1);
00488  goto QUIT;
00489 // BERR2:
00490 // MessageBox(NULL,"Normal - Problem",NULL,MB_OK); goto QUIT;
00491  BERR:
00492  MessageBox(NULL,"Memory Error",NULL,MB_OK);
00493  QUIT:
00494  if(w      != NULL)X__Free(w);
00495  if(Points != NULL)X__Free(Points);
00496  vpp=MainVp; for(i=0;i<Nvert;i++){
00497    if(vpp->status == INEDITOR){
00498      vpp->status=DESELECTED;
00499    }
00500    vpp++;
00501  }
00502  return;
00503 }
00504 

Generated on Sun Apr 27 14:20:10 2014 for OpenFX by  doxygen 1.5.6