DELAUNAY.C

Go to the documentation of this file.

/* --
OpenFX version 2.0 - Modelling, Animation and Rendering Package
Copyright (C) 2000 - 2007 OpenFX Development Team
-- */

/*  DELAUNAY.C  */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <windows.h>

#include "struct.h"
#include "dstruct.h"

#define TRITOP           0
#define TRIFRONT         1
#define TRIRIGHT         2
#define MAXUNIT    4194304L     //1073741824L

#define DLG_CUT                     600
#define DLG_CUT_SPLIT_Y             601
#define DLG_CUT_SPLIT_N             602
#define DLG_CUT_BOTH                603
#define DLG_CUT_WORKPIECE           604


#include "delaunay.h"

#define DESELECTED 0
#define SELECTED   1
#define INEDITOR   3

static HWND      hParent;
static HINSTANCE hThisInstance;


#if __WATCOMC__
int APIENTRY LibMain(HANDLE hDLL, DWORD dwReason, LPVOID lpReserved){
#else
BOOL WINAPI DllMain(HANDLE hDLL, DWORD dwReason, LPVOID lpReserved){
#endif
  switch (dwReason) {
    case DLL_PROCESS_ATTACH: {
      hThisInstance=hDLL;
      if(hDLL == NULL)MessageBox( GetFocus()," NULL instance",NULL, MB_OK);
      break;
    }
    case DLL_PROCESS_DETACH:
      break;
  }
  return TRUE;
}

#if __SC__
#pragma startaddress(DllMain)
#endif

static BOOL CALLBACK BooleanDlgProc(HWND hwnd, UINT msg,
                             WPARAM wparam, LPARAM lparam){
 BOOL err;
 int i;
 static struct bp {BOOL sp,al;} *dp;
 switch( msg ) {
   case WM_INITDIALOG:
     dp=(struct bp *)lparam;
     CentreDlgOnS(hwnd);
     return (TRUE);
   case WM_COMMAND:
      switch(LOWORD(wparam)){
        case IDCANCEL:
          EndDialog(hwnd,FALSE);
          return(TRUE);
        case IDOK:
          EndDialog(hwnd,TRUE);
          return(TRUE);
        default:
          break;
      }
      break;
    default: break;
 }
 return(FALSE);
}

BOOL _Xmodeler
 (HWND parent_window,HWND info_window,X__STRUCTURE *lpevi){
 HCURSOR hSave;
 int type,split;
 struct bp {BOOL sp,al;} d;
 lpEVI=lpevi;
 hParent=parent_window;
 d.sp=TRUE; d.al=TRUE;
 if(DialogBoxParam(hThisInstance,MAKEINTRESOURCE(DLG_CUT),hParent,
             (DLGPROC)BooleanDlgProc,(LPARAM)&d) == FALSE)return FALSE;
 if(d.al)type=0;  else type=1;
 if(d.sp)split=1; else split=0;
 hSave=SetCursor(LoadCursor(NULL,IDC_WAIT));
 if     (ActiveView == TRITOP)IntersectObject(2);
 else if(ActiveView == TRIFRONT)IntersectObject(1);
 else if(ActiveView == TRIRIGHT)IntersectObject(0);
 SetCursor(hSave);
 return TRUE;
}

#define TOL  200000.0  // very important DO not TAMPER
#define TOL1    0.0005 //0.001 // 0.01 changed 4/6/95 for problem with missed point
#define TOL2    0.9995 //0.999 // 0.99 in phase one too near an edge and disregarded
#define TOL3   -1.00000
#define TOL8    1.00000
#define TOL5   -0.001
#define TOL6    1.001
#define TOL9  100.0
#define TOLA  100.0
#define TOLB  10.0
#define TOLC  100.0

static short EdgePierceFace(vector vi,
                            vector p0, vector p1, vector p2,
                            short *code, short *id){
/* code  = 0 point vi is in plane and inside face                */
/*       = 1   "      "      "        on edge id id's the edge   */
/*       = 2   "      "      "        at vertex id id's the edge */
 double ve1,ve2,ve3,vp1,vp2,vp3,q1,q2,q3,det,a,b,ab,
        det1,det2,det3,detmax,dm;
 short k;
  q1=vi[0]-p0[0];  q2=vi[1]-p0[1];  q3=vi[2]-p0[2];
 ve1=p2[0]-p0[0]; ve2=p2[1]-p0[1]; ve3=p2[2]-p0[2];
 vp1=p1[0]-p0[0]; vp2=p1[1]-p0[1]; vp3=p1[2]-p0[2];
 det1=ve1*vp2-vp1*ve2;    /* new way to ensure optimum numerics */
 det2=ve2*vp3-vp2*ve3;
 det3=ve1*vp3-vp1*ve3;
 k=0; detmax=TOL;
 if((dm=fabs(det1)) > detmax){k=1; detmax=dm;}
 if((dm=fabs(det2)) > detmax){k=2; detmax=dm;}
 if((dm=fabs(det3)) > detmax){k=3; detmax=dm;}
 if(k == 0)return 0;
 else if(k == 1){
   a=( vp2*q1-vp1*q2)/det1;
   b=(-ve2*q1+ve1*q2)/det1;
 }
 else if(k == 2){
   a=( vp3*q2-vp2*q3)/det2;
   b=(-ve3*q2+ve2*q3)/det2;
 }
 else if(k == 3){
   a=( vp3*q1-vp1*q3)/det3;
   b=(-ve3*q1+ve1*q3)/det3;
 }
 else return 0;
/* a=1 at p2      b=1 at p1  */
 ab=a+b;
 if(a < TOL5 || a > TOL6 || b < TOL5 || b > TOL6 || ab > TOL6)return 0;
 if(a >= TOL1 && a <= TOL2 && b >= TOL1 && b <= TOL2 && ab <= TOL2){
   *code=0;  /* inside face */
   return 1;
 }
 else if(a < TOL1){ /* along edge 0-1 */
   if(b < TOL1){ /* p0 */
     *code=2; *id=0;
     return 1;
   }
   else if(b > TOL2){ /* p1 */
     *code=2; *id=1;
     return 1;
   }
   else{  /* on edge 0-1 */
     *code=1; *id=0;
     return 1;
   }
 }
 else if(b < TOL1){ /* along edge 0-2 */
   if(a < TOL1){ /* p0 */
     *code=2; *id=0;
     return 1;
   }
   else if(a > TOL2){ /* p2 */
     *code=2; *id=2;
     return 1;
   }
   else{ /* on edge 0-2 */
     *code=1; *id=2;
     return 1;
   }
 }
 else if(ab > TOL2){/* on edge 1-2 */
   *code=1; *id=1;
   return 1;
 }
 return 0;
}

static int O_normalize(vector v){
 double n,nn;
 n= (double)((v[0]*v[0]) + (v[1]*v[1]) + (v[2]*v[2]));
 if(n < 1.e-10)return FAIL;
 nn=sqrt(n);
 n = 1.0 / nn;
 VECSCALE(n,v,v)
 return OK;
}

static short Del_Swap(vector p, long v1l, long v2l, long v3l){
 /* test to see if edge common to triangle pair  v1,v2,v3  and v1,v2,p */
 /* should be swapped so as to keep the triangles as near equilateral  */
 /* as possible.                                                       */
 vertex *v1,*v2,*v3;
 double cosa,cosb,sina,sinb,sinab;
 vector vp3,v13,v23,va,vb,v1p,v2p;
 v1=(MainVp+v1l); v2=(MainVp+v2l); v3=(MainVp+v3l);
 VECSUB((double)v1->xyz,(double)v3->xyz,v13) O_normalize(v13);
 VECSUB((double)v2->xyz,(double)v3->xyz,v23) O_normalize(v23);
 VECSUB((double)v1->xyz,p,v1p) O_normalize(v1p);
 VECSUB((double)v2->xyz,p,v2p) O_normalize(v2p);
 cosa=DOT(v13,v23);
 cosb=DOT(v1p,v2p);
 if(cosa >= 0.0 && cosb >= 0.0)return 0;     /* a < 90 and b < 90 */
 if(cosa <  0.0 && cosb <  0.0)return 1;
 CROSS(v13,v23,va) sina=va[0]*va[0]+va[1]*va[1]+va[2]*va[2];
 if(sina < 1.e-10)return 0; sina = sqrt(sina);
 CROSS(v2p,v1p,vb) sinb=vb[0]*vb[0]+vb[1]*vb[1]+vb[2]*vb[2];
 if(sinb < 1.e-10)return 0; sinb = sqrt(sinb);
 sinab=sina*cosb+sinb*cosa;
 if(sinab < 0.0)return 1;
 return 0;
}

static long BoolAddVertex(vector v, short status){
 vertex *vv;
 if(!UpdateVertexHeap(Nvert+1)){
   MessageBeep(MB_OK);
   return Nvert-1;
 }
 CreateVertex();
 vv=(MainVp+Nvert-1);
 vv->xyz[0]=(long)v[0];
 vv->xyz[1]=(long)v[1];
 vv->xyz[2]=(long)v[2];
 vv->id=0; vv->status=status;
 if(status == DESELECTED){NvertSelect--; NvertDeselect++;}
 return Nvert-1;
}

static FACET *add_small_face(long a, long b, long c,
                             long f, long i, FACET *w,
                             long *nface){
  long nf;
  nf = *nface;
  if((w=(FACET *)X__Realloc(w,(nf+1)*(long)sizeof(FACET)))
       == NULL)return NULL;
  if(!UpdateFaceHeap(Nface+1))return NULL;
  CreateFace(a,b,c);
  CopyFaceProp((MainFp+f),(MainFp+Nface-1));
  w[nf].f=Nface-1;
  (*nface)++;
  return w;
}

static long Del_Edgeid(long v1, long v2, FACET *w, long k,
                       long *fa1, long *fa2, long *e1, long *e2){
 face *f;
 f=(MainFp+w[k].f);
 if((f->V[0] == v1 && f->V[1] == v2) ||
    (f->V[0] == v2 && f->V[1] == v1)){
   *fa1=w[k].fa[1]; *fa2=w[k].fa[2];
   *e1=w[k].ea[1]; *e2=w[k].ea[2];
   return f->V[2];
 }
 else if((f->V[1] == v1 && f->V[2] == v2) ||
         (f->V[1] == v2 && f->V[2] == v1)){
   *fa1=w[k].fa[2]; *fa2=w[k].fa[0];
   *e1=w[k].ea[2]; *e2=w[k].ea[0];
   return f->V[0];
 }
 else if((f->V[2] == v1 && f->V[0] == v2) ||
         (f->V[2] == v2 && f->V[0] == v1)){
   *fa1=w[k].fa[0]; *fa2=w[k].fa[1];
   *e1=w[k].ea[0]; *e2=w[k].ea[1];
   return f->V[1];
 }
 MessageBox(NULL,"Algorithm confused",NULL,MB_OK);
 return -1;
}

static FACET *AddDelaunayTriangulation(vector p2,long *nf, FACET *w){
 short code,id;
 long j,k,l,m,t,n,nn,*Dstack=NULL,Dss=0,fa1,fa2;
 FACET *bf;
 vector pa,pb,pc;
 long vp,w0,w1,w2,vo,vt;
 long f1;
 long e0,e1,e2,ep;
 n= *nf;
 for(j=0;j<n;j++){ /* find which face this intersects    */
   bf=(w+j);
   f1=bf->f;
   w0=(MainFp+f1)->V[0]; w1=(MainFp+f1)->V[1]; w2=(MainFp+f1)->V[2];
   POINT2VECTOR((MainVp+w0)->xyz,pa)
   POINT2VECTOR((MainVp+w1)->xyz,pb)
   POINT2VECTOR((MainVp+w2)->xyz,pc)
   if(EdgePierceFace(p2,pa,pb,pc,&code,&id) /* && code == 0 */){
     vp=BoolAddVertex(p2,DESELECTED);
     if(!UpdateEdgeHeap(Nedge+3))return NULL;
     CreateEdge(vp,w0);  e0=Nedge-1;
     CreateEdge(vp,w1);  e1=Nedge-1;
     CreateEdge(vp,w2);  e2=Nedge-1;
     if((w=add_small_face(vp,w0,w1,f1,j,w,nf))==NULL)return NULL;
     if((w=add_small_face(vp,w2,w0,f1,j,w,nf))==NULL)return NULL;
     (MainFp+f1)->V[0]=vp;       /* modify this face       */
     nn = *nf;
     w[nn-2].fa[0]=nn-1;         /* adjacency in new faces */
     w[nn-2].fa[1]=w[j].fa[0];
     w[nn-2].fa[2]=j;
     w[nn-1].fa[0]=j;
     w[nn-1].fa[1]=w[j].fa[2];
     w[nn-1].fa[2]=nn-2;
     if((k=w[j].fa[0]) >= 0){   /* change adjacency in adjacent */
       for(l=0;l<3;l++)if(w[k].fa[l] == j)w[k].fa[l]=nn-2;
     }
     if((k=w[j].fa[2]) >= 0){
       for(l=0;l<3;l++)if(w[k].fa[l] == j)w[k].fa[l]=nn-1;
     }
     w[j].fa[0]=nn-2;           /*  change adjacency in this face */
     w[j].fa[2]=nn-1;
     w[nn-2].ea[0]=e0;
     w[nn-2].ea[1]=w[j].ea[0];
     w[nn-2].ea[2]=e1;
     w[nn-1].ea[0]=e2;
     w[nn-1].ea[1]=w[j].ea[2];
     w[nn-1].ea[2]=e0;
     w[j].ea[0]=e1;
     w[j].ea[2]=e2;
     if(Dstack == NULL){
       if((Dstack=(long *)X__Malloc(sizeof(long)*(3)))
                 == NULL)return NULL;
     }
     else {
       if((Dstack=(long *)X__Realloc(Dstack,sizeof(long)*(Dss+3)))
                 == NULL)return NULL;
     }
     if(w[j].fa[1] >= 0)   { Dstack[Dss] = j;    Dss++;}
     if(w[nn-2].fa[1] >= 0){ Dstack[Dss] = nn-2; Dss++;}
     if(w[nn-1].fa[1] >= 0){ Dstack[Dss] = nn-1; Dss++;}
     while(Dss > 0){
       Dss--;
       j=Dstack[Dss]; k=w[j].fa[1];
       if(Dss > 0){
         if((Dstack=(long *)X__Realloc(Dstack,sizeof(long)*(Dss))) == NULL){
           return NULL;  /* pop off top of stack */
         }
       }
       if((vo=Del_Edgeid((MainFp+(w[j].f))->V[1],
                         (MainFp+(w[j].f))->V[2],w,k,
                         &fa1,&fa2,&e1,&e2)) < 0)return NULL;
       if(Del_Swap(p2,(MainFp+(w[j].f))->V[1],(MainFp+(w[j].f))->V[2],vo)){
         e0=w[j].ea[1];
         if((m=w[j].fa[2]) >= 0){   /* change adjacency in adjacent */
           for(l=0;l<3;l++)if(w[m].fa[l] == j)w[m].fa[l]=k;
         }
         if(fa1 >= 0){
           for(l=0;l<3;l++)if(w[fa1].fa[l] == k)w[fa1].fa[l]=j;
         }
         (MainEp+e0)->V[0]=vp;
         (MainEp+e0)->V[1]=vo; /* reverse edge */
         vt=(MainFp+(w[j].f))->V[2];
         (MainFp+(w[j].f))->V[2]=vo;   /* update face j */
         w[j].fa[1]=fa1; t=w[j].fa[2]; w[j].fa[2]=k;
         w[j].ea[1]=e1; ep=w[j].ea[2]; w[j].ea[2]=e0;
         (MainFp+(w[k].f))->V[0]=vp;
         (MainFp+(w[k].f))->V[1]=vo;
         (MainFp+(w[k].f))->V[2]=vt; /* face k */
         w[k].fa[0]=j; w[k].fa[1]=fa2; w[k].fa[2]=t;
         w[k].ea[0]=e0; w[k].ea[1]=e2; w[k].ea[2]=ep;
         if((Dstack=(long *)X__Realloc(Dstack,sizeof(long)*(Dss+2))) == NULL){
           return NULL;  /* push both onto stackstack */
         }
         if(w[j].fa[1] >= 0){Dstack[Dss] = j; Dss++;}
         if(w[k].fa[1] >= 0){Dstack[Dss] = k; Dss++;}
       }
     }
     if(Dstack != NULL)X__Free(Dstack);
     return w; /* intersection found */
   }
 }
 return w;
}

static void IntersectObject(int plane){
 /* FACET is the face being worked on           */
 vector *Points=NULL;
 FACET *w=NULL;
 vertex *vpp;
 long w0,w1,w2;
 long i,j,k,l,m;
 vector p0,p1,p2,du,n;
 long n_work_face=0,nPoints=0;
 long xc, d_max_x = -MAXUNIT, d_min_x = MAXUNIT;
 long yc, d_max_y = -MAXUNIT, d_min_y = MAXUNIT;
 long zc, d_max_z = -MAXUNIT, d_min_z = MAXUNIT;
 long radius;
 if(NvertSelect   == 0)return;
 if(NvertDeselect != 0){
   vpp=MainVp; for(i=0;i<Nvert;i++){
     if(vpp->status == DESELECTED){
       vpp->status=INEDITOR;
     }
     vpp++;
   }
 }
 //
 // Make the bounding first triangle
 //
 l=plane;
 vpp=MainVp; for(i=0;i<Nvert;i++){
   if(vpp->status == SELECTED){
     if(l == 2){j=0; k=1; m=2;}
     if(l == 0){j=1; k=2; m=0;}
     if(l == 1){j=2; k=0; m=1;}
     d_min_x=min(d_min_x,vpp->xyz[j]);
     d_max_x=max(d_max_x,vpp->xyz[j]);
     d_min_y=min(d_min_y,vpp->xyz[k]);
     d_max_y=max(d_max_y,vpp->xyz[k]);
     d_min_z=min(d_min_z,vpp->xyz[m]);
     d_max_z=max(d_max_z,vpp->xyz[m]);
   }
   vpp->id=0; vpp++;
 }
 xc = (d_min_x/2 + d_max_x/2);
 yc = (d_min_y/2 + d_max_y/2);
 zc = (d_min_z/2 + d_max_z/2);
 radius = max((d_max_x - d_min_x),(d_max_y - d_min_y));
 n[0]=n[1]=n[2]=0.0;
 if(l == 2){i=0; j=1; k=2; n[2]=1.0;}
 if(l == 0){i=1; j=2; k=0; n[0]=1.0;}
 if(l == 1){i=2; j=0; k=1; n[1]=1.0;}
 p0[i]=(double)xc;
 p0[j]=(double)yc+2.0*(double)radius;
 p0[k]=zc;
 p1[i]=(double)xc-1.732*(double)radius;
 p1[j]=(double)yc-(double)radius;
 p1[k]=zc;
 p2[i]=(double)xc+1.732*(double)radius;
 p2[j]=(double)yc-(double)radius;
 p2[k]=zc;

 BoolAddVertex(p0,DESELECTED); w0=Nvert-1;
 BoolAddVertex(p1,DESELECTED); w1=Nvert-1;
 BoolAddVertex(p2,DESELECTED); w2=Nvert-1;
 if(!UpdateEdgeHeap(Nedge+3))goto QUIT;
 CreateEdge(w0,w1);
 CreateEdge(w1,w2);
 CreateEdge(w2,w0);
 if(!UpdateFaceHeap(Nface+1))goto QUIT;
 CreateFace(w0,w1,w2);
 if((w=(FACET *)X__Malloc((long)sizeof(FACET)))
       == NULL)goto BERR;
 n_work_face=0;
 w[n_work_face].f=Nface-1;
 for(i=0;i<3;i++){
   w[n_work_face].fa[i] = -1;
   w[n_work_face].ea[i] =  i;
 }
 n_work_face++;
 if((Points=(vector *)X__Malloc(NvertSelect*(long)sizeof(vector)))
      == NULL)goto BERR;
 nPoints=0;
 for(vpp=MainVp,i=0;i<Nvert;i++){
   if(vpp->status == SELECTED){
     POINT2VECTOR(vpp->xyz,Points[nPoints])
     nPoints++;
   }
   vpp++;
 }
 for(i=0;i<nPoints;i++){
   double mu;
   VECSUB(Points[i],p0,du)
   mu=DOT(du,n);
   VECSCALE(mu,n,du)
   VECSUM(Points[i],du,p2)
   if((w=AddDelaunayTriangulation(p2,&n_work_face,w)) == NULL)goto BERR;
 }
 (MainVp+w0)->status=SELECTED;
 (MainVp+w1)->status=SELECTED;
 (MainVp+w2)->status=SELECTED;
 NvertSelect   += 3;
 NvertDeselect -= 3;
 EraseVertex(1);
 goto QUIT;
// BERR2:
// MessageBox(NULL,"Normal - Problem",NULL,MB_OK); goto QUIT;
 BERR:
 MessageBox(NULL,"Memory Error",NULL,MB_OK);
 QUIT:
 if(w      != NULL)X__Free(w);
 if(Points != NULL)X__Free(Points);
 vpp=MainVp; for(i=0;i<Nvert;i++){
   if(vpp->status == INEDITOR){
     vpp->status=DESELECTED;
   }
   vpp++;
 }
 return;
}