/* * glmesh2.c: 2-D triangle mesh generation using OpenGL. * * This program reads a file containing a 2-D triangle mesh, allows * an interactive user to update the mesh with additions, deletions, * and alterations, and writes a file containing the revised mesh. * * R. Renka * renka@cs.unt.edu * * 06/22/2009 * * The required computing platform is a workstation running OpenGL * and the OpenGL Utility Toolkit GLUT or one of its alternatives. * In the case of Apple OS X, both OpenGL and GLUT are installed as * part of Xcode. OpenGL is included in most Linux distributions, * and is installed with Microsoft Visual Studio, but it may be * necessary to download GLUT: * * http://www.opengl.org/resources/libraries/glut/ * * Three files are needed: a static link library (libglut.a or * glut32.lib in the case of Windows), a dynamic link library * (libglut.so.3.8.0, along with soft links libglut.so.3 and * libglut.so or, for Windows, glut.dll), and a header file glut.h. * * Compile and link command for Linux/X11/gcc: * * gcc -O3 glmesh2.c -o glmesh2 -lglut -lGLU -lGL -lXmu -lXi \ * -lXext -lX11 -lm * * Compile and link command for Apple OS X/Cocoa: * * gcc -O3 glmesh2.c -o glmesh2 --ansi -framework GLUT * -framework OpenGL -framework Cocoa * * Compile and link command for Windows Visual C++ (executed in a * command window obtained from Visual Studio so that the path * environment variable is appropriately set): * * cl glmesh2.c * * Program execution command: * * glmesh2 input_file output_file * * A sample input file is provided with this source code file. Its * format is described in strings filemsg1 and filemsg2, which appear * on the screen if the program is executed with fewer than two file * names specified on the command line. Once the program is started, * a menu of options is displayed by a right mouse button press in * the window. * * Note that the ordering of header files is not arbitrary. In * particular, any subset of the following must appear in the order * specified: windows.h, stdlib.h, glew.h, glut.h. * */ #include #include #include #include #ifdef __APPLE__ #include #else #include #endif /* * Structures: * * Rectangles are used in a dialog window. */ struct rect { int x0, y0; /* Window coordinates of lower left corner */ int w, h; /* Width and height in window coordinates */ }; /* * A queue of boundary edges Qe and a priority queue of triangles Qt, * with pointers Qtp, are used by the Delaunay refinement function * refineRd. The vertex indices are stored in order to determine if * a queued object is still in the triangle mesh when it is dequeued. */ struct boundary_edge { /* Boundary edge in queue Qe */ int index_e; /* Edge index */ int index_v1; /* First endpoint vertex index */ int index_v2; /* Second endpoint vertex index */ }; struct triangle { /* Triangle in priority queue Qt */ int index_qtp; /* Index of the pointer in Qtp */ int index_t; /* Triangle index */ int index_v1; /* First vertex index */ int index_v2; /* Second vertex index */ int index_v3; /* Third vertex index */ int imin; /* Local index of vertex opposite shortest edge */ double emin; /* Squared length of shortest edge */ }; /* * Static array dimension: */ #define LSTRING 80 /* Length of title string and dialog string */ /* * Enumeration constant and arrays of strings for printing. */ enum Mode {DEFAULT, SELECT_CIRCLE, PAN, CONSTRUCT_R, ALTER_R, ALTER_VIEW_VOLUME, ALTER_DOMAIN, SWAP_EDGES, PLACE_TITLE, EDIT_UV}; static char *modeName[] = {"DEFAULT", "SELECT_CIRCLE", "PAN", "CONSTRUCT_R", "ALTER_R", "ALTER_VIEW_VOLUME", "ALTER_DOMAIN", "SWAP_EDGES", "PLACE_TITLE", "EDIT_UV"}; static char *flagName[] = {"FALSE", "TRUE"}; /* * Global variables defining the triangle mesh: * * A triangle mesh is a set of triangles in which the intersection * of two triangles, if not empty, is a vertex of both or an edge * of both, and whose union is a connected region with one or more * simple closed non-intersecting boundary curves. Boundary curves * are oriented so that triangles are on the left as the sequence * of boundary vertices is traversed. The domain geometry can be * defined by designating a subset of the boundary vertices as being * non-removable. * * A Delaunay triangulation is a triangle mesh in which the triangle * circumcircles contain no vertices in their interiors. A * constrained Delaunay triangulation (CDT) is a triangle mesh in * which, if the circumcircle of a triangle has a vertex in its * interior, that vertex is separated from (every point of) the * triangle interior by a boundary edge. Refer to Functions * optimizeRt, refineRd, swap, and swptst. * * The following variables characterize the mesh. * * nvmax = Maximum number of vertices for which storage has been * allocated. This value is increased, and additional * storage is allocated when necessary as the mesh is * updated with additional vertices. * * nc >= 1: Number of boundary curves. * nb >= 3*nc: Number of boundary vertices (boundary edges). * nv >= nb: Number of vertices. * nn <= nb: Number of non-removable vertices. * nt = 2*nv - nb + 2*nc - 4 <= 2*nv - 5: Number of triangles. * ne = 3*nv - nb + 3*nc - 6 <= 3*nv - 6: Number of edges. * * vtxy = Array dimensioned nvmax by 2 containing the Cartesian * coordinates of the vertices. * * uv = Array dimensioned nvmax by 2 containing pairs of function * values (u,v), such as velocity vectors, at the vertices. * The values default to zeros, but may be altered. * * ltri = Array dimensioned 2*nvmax by 9 containing a triangle list * in the first nt rows: CCW-ordered vertex indices (0 to * nv-1), neighboring triangle indices (-1 to nt-1), and * opposite edge indices (0 to ne-1), where nonexistent * neighboring triangles are indexed by -1. The ordering of * the neighboring triangle indices and edge indices is * defined by that of the vertices: edge ltri[kt][i+6] is * opposite vertex ltri[kt][i] and is shared by triangles kt * and ltri[kt][i+3] for i = 0, 1, and 2. * * lvtx = Array of length nvmax containing a vertex list in the * first nv positions: index of a triangle that contains * vertex kv in position kv for kv = 0 to nv-1. If kv indexes * a boundary vertex, then triangle lvtx[kv] also contains the * next vertex in the boundary curve, where boundary curves * are oriented so that triangles are on the left. * * ledg = Array of length 3*nvmax containing an edge list in the * first ne positions: index of a triangle (the one with * larger index if there are two) that contains edge ke in * position ke for ke = 0 to ne-1. * * lnrv = Array of length nvmax containing the indices (for vtxy) of * the non-removable vertices in the first nn positions. * * lcrv = Array of length nvmax/3 containing the index of one * boundary vertex in each boundary curve in the first nc * positions. * * lste = Temporary storage array of length 3*nvmax used by several * functions to store lists of edges, triangles, or vertices, * and used by Function display to store triangle color * indices. * * listn = Array of length nvmax used to store Boolean flags * marking vertices as non-removable. This array must be * compatible with the contents of lnrv. */ static GLint nvmax, nc, nb, nv, nn = 0, nt, ne; static GLdouble (*vtxy)[2], (*uv)[2]; static GLint (*ltri)[9], *lvtx, *ledg, *lnrv, *lcrv, *lste; static GLboolean *listn; /* * Global variables associated with the currently selected region R: * * The user may select a polygonal region R to be refined, optimized, * coarsened, or removed from the mesh. The polygon boundary is a * simple closed unoriented curve defined by a sequence of three or * more vertices (selected by mouse button presses). The following * variables are associated with R. * * mode = Current mode (described below). Values CONSTRUCT_R and * ALTER_R are related to R. * newR = Flag with value TRUE iff Function initR needs to be called * to store flags listR before the mesh can be altered or * correctly displayed: set to TRUE when R is created and * when R is altered by adding, deleting, or moving a vertex * (Function mouse); set to FALSE by Function initR; and * tested by Functions coarsenRv, convertRv, fillRt, moveRn, * optimizeRt, optimizeRv, refineRd, refineRe, refineRt, * removeRt, and trprint. * indxR = Index (for R, in the range 0 to np-1) of a vertex * currently being dragged by the mouse. * listR = Array of length nvmax used to store Boolean flags * marking vertices as contained in R. An edge is contained * in R if and only if its endpoint vertices are in R, and * a triangle is in R if and only if its vertices are in R. * np = Number of vertices in the boundary of the currently selected * polygonal region R. * npmax = Maximum value of np for which it is not necessary to * reallocate storage for R. * R = Array dimensioned npmax by 2 containing the Cartesian coordi- * nates of the sequence of vertices defining the boundary of R. */ static GLboolean newR = 1; static GLint indxR; static GLboolean *listR; static GLint np = 0; static GLint npmax = 100; static GLdouble (*R)[2]; /* * Additional global variables: * * antialias = Flag with value TRUE iff edges are to be drawn with * antialiasing. * aspect = Flag with value TRUE iff the viewport is required to * have the same aspect ratio as the viewing volume. * bdry_layer = Flag with value TRUE iff the smoothing function * optimizeRv is to create a boundary layer by moving * vertices adjacent to boundary vertices closer to * the boundary. * color_b = Color of background. * color_e = Color of triangle mesh edges. * color_fill = Flag specifying the method used by Function display * to render triangles: * color_fill = 0 if triangles are rendered in outline * mode by drawing only the edges * (using color_e). * color_fill = 1 if all triangles are to be color- * filled using four alternating colors * (a 4-coloring of the mesh using array * colors) in addition to edge drawing. * color_fill = 2 if only the triangles contained in R * are to be color-filled as above. * color_fill = 3 if triangles are to be color-coded by * quality using function trqual1. * color_t = Color of text strings. * colors = Array dimensioned 4 by 3 containing the triangle colors * used when color_fill is 1 or 2. * cr,cx,cy = Radius and center of a circle selected by Functions * getCircle and motion. * dlg_exit = Flag with value TRUE iff the program is to be * terminated when the user selects the 'Yes' button in * the dialog box. * dlg_r0 = Rectangle representing 'No' button. * dlg_r1 = Rectangle representing 'Yes' button. * dlg_string = Text string displayed in dialog box window. * dlg_yes = Flag with value TRUE iff the user selects the 'Yes' * button in the dialog box. * dx,dy = Width and height of viewing volume. * fname = File name for output. * ind_font = Font number for indices. * line_width = Line width in pixels for drawing edges. * min_angle = Minimum angle (in degrees) in a triangle for which the * triangle is acceptable: criterion used by Functions * badTriangle and refineRd. * mode = Current mode, altered by a keypress or menu selection, and * tested by Function motion: * mode = SELECT_CIRCLE if the user is selecting a circle. * Refer to Function moveRn. * mode = PAN if the user is selecting a new center (xc,yc). * mode = CONSTRUCT_R if the user is selecting a new polygonal * region R. * mode = ALTER_R if the user is altering polygonal region R. * mode = ALTER_VIEW_VOLUME if the user is selecting a new * viewing volume by dragging a rectangle. * mode = ALTER_DOMAIN if the user is altering the domain by * dragging non-removable vertices. * mode = SWAP_EDGES if the user is selecting triangle edges * in convex quadrilaterals for edge swaps. * mode = PLACE_TITLE if the user is selecting a title * location. * neq = Number of boundary edges in Qe. * neqmax = Maximum value of neq for which storage has been * allocated. * nnbv = Number of neighbors of a non-removable vertex that is * currently being dragged with mode = ALTER_DOMAIN. * nref = Number of refinements (by refineRd, refineRe, or refineRt) * that have not been reversed by Function unrefine, or * rendered non-reversible by Function coarsenRv or removeRt. * nrefmax = Maximum value of nref for which it is not necessary to * reallocate storage for rstack. * nsp = Number of Steiner points (vertices added to the triangle * mesh) in the previous call to Function refineRd. * nspmax = Maximum number of Steiner points that may be used by each * call to Function refineRd. This maximum will be reached * if min_angle is too large. * ntq = Number of triangles in Qt. * ntqmax = Maximum value of ntq for which storage has been * allocated. * offcenter = Flag with value TRUE iff the Delaunay refinement * method, Function refineRd, is to use off-centers, * rather than circumcenters, as Steiner points. * offscale = Scale factor for the minimum edge length, used to * compute the distance from the edge to an off-center: * 0.96*(1+c)/(2*s) = 0.48*sqrt((1+c)/(1-c)), where * c and s are the cosine and sin of min_angle. The * initial value is computed by Function storeV. * plot_box = Flag with value TRUE iff an axis-aligned bounding box * with labeled corner coordinates is to be drawn. * plot_e = Flag with value TRUE iff edge indices are to be displayed * (at edge centers). * plot_nrv = Flag with value TRUE iff non-removable vertex indices * are to be displayed (with color distinct from that of * vertex indices). * plot_t = Flag with value TRUE iff triangle indices are to be * displayed (near barycenters). * plot_v = Flag with value TRUE iff vertex indices are to be * displayed. * Qe = Array of length neqmax containing a queue of struct boundary_ * edge objects used by Function refineRd to store boundary * edges for splitting. * qefront,qeback = Indices of the front and back of the queue in * circular array Qe. * Qt = Array of length ntqmax containing a priority queue of * struct triangle objects used by Function refineRd to * prioritize triangles for splitting. * Qtp = Array of length ntqmax+1 containing a binary heap of * pointers to triangles in Qt. * rstack = Array of length nrefmax containing a sequence of vertex * counts (nv values) associated with refinements by * Functions refineRd, refineRe, and refineRt, and unre- * finements by Function unrefine. * scamax = Squared cosine of min_angle. The initial value is * computed in Function storeV. * sx,sy = Scale factors in the mapping from window coordinates * to object coordinates -- computed in reshape. * title = Optional character string displayed in the primary window. * title_font = Font number for the title string. * title_pos = Raster position of (lower left corner of) title. * window = Array of window identifiers. * xc,yc = Object coordinates of the point that maps to the center * of the viewport: initialized to the center of the * bounding box, but may be altered to pan around (when * zoomed in, for example). * xmin,ymin = Coordinates of the lower left corner of the bounding * box containing the vertices. * xmax,ymax = Coordinates of the upper right corner of the bounding * box containing the vertices. * xv1,yv1 = Object coordinates of the lower left or upper right * corner of the viewing volume, or of an endpoint of a * line segment during a drag operation: used by Function * motion. * xv2,yv2 = Object coordinates of the upper right or lower left * corner of the viewing volume, or of the cursor position * during a drag operation: used by Function motion. * xv3,yv3 = Object coordinates of an endpoint of a line segment * during a drag operation: used by Function motion. */ static GLboolean antialias = GL_FALSE; static GLboolean aspect = GL_TRUE; static GLboolean bdry_layer = GL_FALSE; static GLfloat color_b[3] = {0.0, 0.0, 0.0}; static GLfloat color_e[3] = {1.0, 1.0, 1.0}; static GLint color_fill = 0; static GLfloat color_t[3] = {1.0, 1.0, 1.0}; static GLdouble cr = 0.0, cx, cy; static GLboolean dlg_exit = GL_FALSE; struct rect dlg_r0 = {240, 25, 90, 35}; struct rect dlg_r1 = {70, 25, 90, 35}; static char dlg_string[LSTRING] = ""; static GLboolean dlg_yes = GL_FALSE; static GLdouble dx, dy; static char *fname; static unsigned int ind_font = 2; static GLfloat line_width = 2.0; static int min_angle = 20; static enum Mode mode = DEFAULT; static int neq = 0; static int neqmax; static GLint nnbv = 0; static GLint nref = 0; static GLint nrefmax = 10; static int nsp = 0; static int nspmax = 300; static int ntq = 0; static int ntqmax; static GLboolean offcenter = GL_TRUE; static double offscale; static GLboolean plot_box = GL_FALSE; static GLboolean plot_e = GL_FALSE; static GLboolean plot_nrv = GL_TRUE; static GLboolean plot_t = GL_FALSE; static GLboolean plot_v = GL_FALSE; static struct boundary_edge *Qe = NULL; static int qefront, qeback; static struct triangle *Qt = NULL; static struct triangle **Qtp = NULL; static GLint *rstack; static double scamax; static GLdouble sx, sy; static char title[LSTRING] = ""; static unsigned int title_font = 3; static GLfloat title_pos[3] = {0., 0., 0.}; static int window[2]; static GLdouble xc, yc; static GLdouble xmin, ymin, xmax, ymax; static GLdouble xv1, yv1, xv2, yv2, xv3, yv3; /* * Color Table for rendering triangles: */ static GLfloat colors[4][3] = { {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}, {1.0, 1.0, 0.0} }; /* * Function descriptions: * * adjacent: Test for existence of an edge between a pair of * vertices. * * badTriangle: Test a triangle for a small angle (defined by * global variables min_angle and scamax). * * bbox: Compute and store the corner coordinates of the bounding * box: [xmin,xmax] X [ymin,ymax]. * * bcurve: Return the CCW-ordered sequence of boundary vertex * indices in a boundary curve. * * buildQt: Construct a priority queue Qt containing all the bad * triangles in R, as defined by Function badTriangle, * with priority given to shortest edge length. * * closeR: Close a user-selected polygon R by implicitly adding * an edge from the last vertex to the first. * * coarsenRv: Coarsen R by removing about 70% of the vertices. * * convertRv: Convert boundary vertices in R to removable or non- * removable by deleting their indices from lnrv or * inserting their indices into lnrv, and adjusting * the flags in listn accordingly. * * deleteQt: Return the first triangle in priority queue Qt, and * remove it from the queue. * * deltri: Delete a triangle from the mesh. * * delvtx: Delete a vertex from the mesh. * * display: Draw triangle mesh edges with optional triangle fill, * index annotation, title, etc. * * displayPosition: Display the coordinates of a point, such as the * current mouse position, in the upper left corner * of the viewport. * * displayString: Write a string to the screen. * * dlgDisplay: Dialog window display callback. * * dlgKey: Keyboard callback associated with the dialog window. * * dlgMouse: Callback triggered by a mouse button press or release * with the mouse position in the dialog window. * * dlgReshape: Reshape function for the dialog window (window 1). * * dragCorner: Mouse callback used to alter the domain by moving * a non-removable vertex. * * drawCircle: Draw a circle of radius cr centered at (cx,cy). * * drawRect: Draw the outline of an axis-aligned rectangle. * * encroached: Test a boundary edge for encroachment: the vertex * opposite the edge inside its diametral circle. * * fillRt: Fill R with triangles where R contains three adjacent * boundary vertices forming a reentrant angle (corner) at * a removable vertex. * * filterDown: Filter down an insertion point in a binary heap in * accordance with the heap-order property: used by * buildQt and deleteQt. * * getCenter: Mouse callback used to move the center (xc,yc). * * getCircle: Mouse callback used to select a circle center (cx,cy) * and radius cr. * * getCorners: Replace the set of non-removable vertices lnrv with * the set of boundary vertices that form corners (do * not lie on the lines defined by their predecessors * and successors), and adjust the flags in listn. * * getLin: Read a line from stdin. * * getPolygon: Mouse callback used to construct R. * * getTitlePos: Mouse callback used to set the title location. * * getVertex: Mouse callback used to select a vertex in uv edit * mode. * * getViewVolume: Mouse callback used to change the viewing volume * by dragging a rectangle. * * init0: Initialization for window 0. * * initR: Initialization for user-selected polygon R: store flags * listR. * * inputData: Allocate storage and read a data set containing a tri- * angle mesh: nv, vtxy, nt, and first three columns of * ltri. * * insertQt: Insert a triangle into priority queue Qt. * * insertv: Partition a triangle into three subtriangles by insert- * ing a new vertex in its interior. * * inside: Locate a point (xp,yp) relative to the polygon R, * returning a nonzero value iff (xp,yp) is contained in R. * * intsec: Test for an intersection between a pair of line * segments in the plane. * * key: Keyboard callback for window 0 (calls menu). * * makeMenu: Create a menu for window 0. * * menu: Respond to menu selections. * * motion: Callback for mouse motion with a button pressed * * mouse: Default mouse callback: used to alter R. * * moveRn: Project all non-removable vertices in R (that can be * moved without destroying the mesh) onto the current * circle (if any). * * offCenter: Compute an off-center of a triangle: a point on the * perpendicular bisector of the shortest edge, between * the midpoint of the edge and the circumcenter of the * triangle with distance from the edge just small enough * that the new triangle formed by the edge and the off- * center will not be selected for splitting by Function * badTriangle. * * optimizeRt: Optimize R with Delaunay edge swaps in convex * quadrilaterals consisting of pairs of adjacent * triangles in R. * * optimizeRv: Optimize R by applying a single iteration of a * smoothing operation in which each removable vertex * is moved (if possible) to a weighted sum of the * vertices in its star (the triangles containing the * vertex). * * passiveMotion: Callback triggered by mouse motion with no mouse * buttons pressed: used to display the object * coordinates of the mouse position. * * plDistance: Test for a point-line distance within a tolerance. * * reallocQe: Reallocate storage for Qe, and copy the contents of * the old array into the new array with an ordering * that reflects the circular array storage scheme. * * refineRd: Refine R by Delaunay refinement. Bad triangles are * split by inserting their circumcenters or off-centers, * and encroached boundary edges are split by inserting * their midpoints (or alternative split points). Edges * are then swapped as required by the empty circumcircle * criterion (Function swptst). If a constrained Delaunay * triangulation is input, a Delaunay triangulation will * be output. * * refineRe: Refine R by adding a new vertex at the midpoint of each * edge in R (using split and swap), partitioning each * triangle into four equal-area subtriangles. * * refineRt: Refine R by adding triangle barycenters, partitioning * each triangle into three equal-area subtriangles. * * removeRt: Remove the triangles in R, creating new boundary curves * if necessary. * * reshape: Reshape callback (choose viewport and viewing volume) * for window 0. * * specialKey: Move center (xc,yc) in response to arrow keys, zoom * in or out on Page Up or Page Down keys, or restore * default viewing volume if Home key is pressed. * * split: Add a new vertex in the interior of an edge, splitting the * edge into a pair of edges, and replacing each triangle * containing the edge by two triangles. * * splitEdges: Split all encroached boundary edges in R by inserting * their midpoints (or alternative split points) into * the triangle mesh: used by Function refineRd. * * storeV: Compute and store global parameter values that depend on * the input data: dx, dy, xc, xmax, xmin, yc, ymax, and * ymin, and allocate storage for R and rstack. * * swap: Modify a triangle mesh by swapping diagonal edges in a con- * vex quadrilateral defined by a pair of adjacent triangles. * * swapEdges: Mouse callback allowing the user to swap edges. * * swptst: Delaunay circumcircle test. * * trealloc: Increase nvmax, and re-allocate storage for the * triangle mesh. * * trfind: Locate a point relative to the triangle mesh. * * trmtst: Triangulation data structure integrity test. * * trprint: Print triangle mesh data structure and statistics on * stdout. * * trqual1: Compute a measure of triangle quality: triangle area * divided by mean side length, normalized to [0,1]. * * trqual2: Compute a measure of triangle quality: twice the ratio * of the inradius to the circumradius. * * trwrite: Output the triangle mesh to a file. * * unrefine: Unrefine the mesh by reversing the previous refinement * by refineRd, refineRe, or refineRt, if possible. * * unswap: Reverse the effect of a call to Function swap: swap * diagonal edges in a convex quadrilateral defined by a * pair of adjacent triangles, with the choice of new * triangle indices reversed from that of Function swap. * * update: Update a Delaunay triangulation in the vicinity of a * newly inserted vertex by applying a sequence of swap * tests and the appropriate swaps to the edges opposite * the new vertex. The queues Qe and Qt are also updated. * * xlate: Translate glut window coordinates (wx,wy), with depth wz, * to object coordinates (x,y,z) by inverting the mappings * associated with the modelview, projection, and viewport * matrices. * * main: Input data, create windows, create a menu, register call- * back functions, and enter the event loop under control of * glut. */ /* * Function prototypes */ static int adjacent(int kv1, int kv2); static int badTriangle(int kt, int *imin, double *emin, double *scamx); static void bbox(void); static unsigned char bcurve(int kc, int *nbv, int indbv[]); static void buildQt(void); static void closeR(void); static void coarsenRv(void); static void convertRv(unsigned char removable); static void deleteQt(struct triangle *ptri); static unsigned char deltri(int kt0); static unsigned char delvtx(int kv0); static void display(void); static void displayPosition(GLdouble x, GLdouble y); static int displayString(GLfloat pos[3], char *string, unsigned int fontno, GLfloat scolor[3]); static void dlgDisplay(void); static void dlgKey(unsigned char key, int x, int y); static void dlgMouse(int button, int state, int x, int y); static void dlgReshape(int w, int h); static void dragCorner(int button, int state, int x, int y); static void drawCircle(void); static void drawRect(GLuint x0, GLuint y0, GLuint w, GLuint h, GLfloat brcolor[3], GLfloat tlcolor[3]); static int encroached(int ke, int *iv1, int *iv2); static void fillRt(void); static void filterDown(int ip); static void getCenter(int button, int state, int x, int y); static void getCircle(int button, int state, int x, int y); static void getCorners(void); static int getLin(char *s, int len); static void getPolygon(int button, int state, int x, int y); static void getTitlePos(int button, int state, int x, int y); static void getVertex(int button, int state, int x, int y); static void getViewVolume(int button, int state, int x, int y); static void init0(void); static void initR(void); static void inputData(int argc, char *argv[]); static void insertQt(int kt, int imin, double emin); static void insertv(int kt0, double b1, double b2, double b3); static unsigned char inside(double xp, double yp); static unsigned char intsec(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4); static void key(unsigned char key, int x, int y); static void makeMenu(void); static void menu(int item); static void motion(int x, int y); static void mouse(int button, int state, int x, int y); static void moveRn(void); static void offCenter(int kt, int imin, double *cx, double *cy); static void optimizeRt(void); static void optimizeRv(void); static void passiveMotion(int x, int y); static unsigned char plDistance(double xp, double yp, double x1, double y1, double x2, double y2, double tols); static void reallocQe(void); static void refineRd(void); static void refineRe(void); static void refineRt(void); static void removeRt(void); static void reshape(int w, int h); static void specialKey(int key, int x, int y); static void split(int ke, double b1, double b2); static void splitEdges(unsigned char tritest); static void storeV(void); static int swap(int ke); static void swapEdges(int button, int state, int x, int y); static unsigned char swptst(int in1, int in2, int io1, int io2); static void trealloc(int nvnew); static int trfind(int kt0, double xp, double yp, int *kep, int *ktp, double *b0p, double *b1p, double *b2p); static void trmtst(int nc, int nb, int nv, int nn, int nt, int ne, double vtxy[][2], int ltri[][9], int lvtx[], int ledg[], int lnrv[], int lcrv[], unsigned char listn[]); static void trprint(int nc, int nb, int nv, int nn, int nt, int ne, double vtxy[][2], double uv[][2], int ltri[][9], int lvtx[], int ledg[], int lnrv[], int lcrv[]); static double trqual1(double x1, double y1, double x2, double y2, double x3, double y3); static double trqual2(double x1, double y1, double x2, double y2, double x3, double y3); static void trwrite(char *fname, int nc, int nb, int nv, int nt, int ne, double vtxy[][2], double uv[][2], int ltri[][9]); static void unrefine(void); static void unswap(int ke); static void update(int kv, unsigned char tritest, int *ns, int inde[]); static void xlate(GLint wx, GLint wy, GLdouble wz, GLdouble *x, GLdouble *y, GLdouble *z); /* *-------------------------------------------------------------------- * * adjacent: Test for existence of a triangulation edge between a * pair of vertices. * * On input: kv1 and kv2 are vertex indices in the range 0 to nv-1. * * On output: the return value is 1 if the vertices are adjacent, * or 0 otherwise. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static int adjacent(int kv1, int kv2) { int i, kt, kt0, kts; /* * Loop on triangles kt containing vertex kv1, beginning with kt0. */ kt0 = lvtx[kv1]; kt = kt0; do { if (ltri[kt][0] == kv1) i = 1; else if (ltri[kt][1] == kv1) i = 2; else i = 0; if (ltri[kt][i] == kv2) return 1; kts = kt; kt = ltri[kt][i+3]; } while (kt >= 0 && kt != kt0); /* * Test the last neighbor if kv1 is a boundary vertex. */ if (kt < 0) { if (ltri[kts][(i+1)%3] == kv2) return 1; } /* * kv2 is not a neighbor of kv1. */ return 0; } /* End of adjacent */ /* *-------------------------------------------------------------------- * * badTriangle: Test a triangle for a small angle (defined by global * variables min_angle and scamax). * * This function is used by Function refineRd to select triangles * for splitting. * * On input: kt is the index (0 to nt-1) of the triangle to be * tested. * * On output: imin is the local index (ltri column index) of the * vertex opposite the shortest edge in triangle kt, * emin is the squared length of the shortest edge in * triangle kt, scamx is the squared cosine of the * smallest angle in triangle kt, and the return value * is 1 if the smallest angle in triangle kt is less * than min_angle (and does not separate two boundary * edges), or 0 otherwise. * * Note that this function could be modified to test additional * criteria, such as triangle size or error estimates associated * with triangles, returning 1 if any of the conditions for splitting * are met. * * Note also that small angles defined by non-removable vertices * cannot be eliminated. * * Global variables required: ltri, scamax, vtxy * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static int badTriangle(int kt, int *imin, double *emin, double *scamx) { double dsmin, ds1, ds2, ds3, dx1, dx2, dx3, dy1, dy2, dy3, sca1, sca2, sca3, scam, x1, x2, x3, y1, y2, y3; int i, kv1, kv2, kv3; kv1 = ltri[kt][0]; kv2 = ltri[kt][1]; kv3 = ltri[kt][2]; x1 = vtxy[kv1][0]; y1 = vtxy[kv1][1]; x2 = vtxy[kv2][0]; y2 = vtxy[kv2][1]; x3 = vtxy[kv3][0]; y3 = vtxy[kv3][1]; dx1 = x3-x2; dy1 = y3-y2; dx2 = x1-x3; dy2 = y1-y3; dx3 = x2-x1; dy3 = y2-y1; ds1 = dx1*dx1 + dy1*dy1; ds2 = dx2*dx2 + dy2*dy2; ds3 = dx3*dx3 + dy3*dy3; i = 0; dsmin = ds1; if (ds2 < dsmin) { i = 1; dsmin = ds2; } if (ds3 < dsmin) { i = 2; dsmin = ds3; } *imin = i; *emin = dsmin; /* * The smallest angle is less than min_angle iff its squared * cosine (scamx) is greater than that of min_angle (scamax). */ if (dsmin == 0.0) return 1; sca1 = dx3*dx2 + dy3*dy2; sca1 = sca1*sca1/(ds3*ds2); sca2 = dx1*dx3 + dy1*dy3; sca2 = sca2*sca2/(ds1*ds3); sca3 = dx2*dx1 + dy2*dy1; sca3 = sca3*sca3/(ds2*ds1); i = 0; scam = sca1; if (sca2 > scam) { i = 1; scam = sca2; } if (sca3 > scam) { i = 2; scam = sca3; } if (scam > 1.0) scam = 1.0; *scamx = scam; return ((ltri[kt][(i+1)%3+3] >= 0 || ltri[kt][(i+2)%3+3] >= 0) && scam > scamax); } /* End of badTriangle */ /* *-------------------------------------------------------------------- * * bbox: Compute and store the corner coordinates defining the * bounding box. * * Global variables required: nv, vtxy, xmin, xmax, ymin, ymax * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void bbox(void) { int i; xmin = vtxy[0][0]; xmax = xmin; ymin = vtxy[0][1]; ymax = ymin; for (i = 1; i < nv; i++) { if (vtxy[i][0] < xmin) xmin = vtxy[i][0]; if (vtxy[i][0] > xmax) xmax = vtxy[i][0]; if (vtxy[i][1] < ymin) ymin = vtxy[i][1]; if (vtxy[i][1] > ymax) ymax = vtxy[i][1]; } return; } /* End of bbox */ /* *-------------------------------------------------------------------- * * bcurve: Return the CCW-ordered sequence of boundary vertex * indices in a boundary curve. * * On input: kc is a curve index in the range 0 to nc-1. * * On output: The return value is 0 if kc is invalid. Otherwise, * the return value is 1, nbv is the number of vertices * in the curve, and indbv contains the sequence of * indices in the first nbv positions. * * Global variables required: lcrv, ltri, lvtx, nc * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static unsigned char bcurve(int kc, int *nbv, int indbv[]) { int i, kt, kv, kv0, lbv; /* * Test for invalid input. */ if (kc < 0 || kc >= nc) return 0; /* * Loop on boundary vertices kv starting with kv0. */ lbv = 0; kv0 = lcrv[kc]; kv = kv0; do { indbv[lbv] = kv; lbv++; kt = lvtx[kv]; if (ltri[kt][0] == kv) i = 1; else if (ltri[kt][1] == kv) i = 2; else i = 0; kv = ltri[kt][i]; } while (kv != kv0); *nbv = lbv; return 1; } /* End of bcurve */ /* *-------------------------------------------------------------------- * * buildQt: Allocate storage and construct a priority queue Qt * containing all the bad triangles in R, as defined by * Function badTriangle, with priority given to shortest * edge length. * * The priority queue is implemented as a binary heap: a complete * binary tree of pointers stored in array Qtp positions 1 to ntq * so that the children of Qtp[i] are in Qtp[2*i] and Qtp[2*i+1], * and the parent is in Qtp[i/2], with the heap-order property: * every node has key value (emin value) less than or equal to those * of its children. * * Related function deleteQt removes the top-priority triangle at a * cost of O(log(ntq)) operations to restore heap-order property, * and insertQt inserts a new triangle in constant expected time * (but worst-case logarithmic time). Rather than using insertQt, * this function inserts the triangles in the order in which they are * encountered, and then employs a linear-time algorithm to restore * the heap-order property. * * Global variables required: ltri, ntq, ntqmax, Qt, Qtp, scamax, * vxtx * * glmesh2 functions called: badTriangle, filterDown * *-------------------------------------------------------------------- */ static void buildQt(void) { double emin, scamx; int imin, ip, kt; ntqmax = nt; ntq = 0; /* * Allocate storage for Qt and Qtp: enough for all triangles in * the mesh. */ Qt = malloc(ntqmax*sizeof *Qt); Qtp = malloc((ntqmax+1)*sizeof *Qtp); /* * Loop on triangles kt in R, storing the bad ones in Qtp without * enforcing the heap-order property. */ for (kt = 0; kt < nt; kt++) { if (!listR[ltri[kt][0]] || !listR[ltri[kt][1]] || !listR[ltri[kt][2]]) continue; if (badTriangle(kt, &imin, &emin, &scamx)) { Qt[ntq].index_qtp = ntq+1; Qt[ntq].index_t = kt; Qt[ntq].index_v1 = ltri[kt][0]; Qt[ntq].index_v2 = ltri[kt][1]; Qt[ntq].index_v3 = ltri[kt][2]; Qt[ntq].imin = imin; Qt[ntq].emin = emin; ntq++; Qtp[ntq] = &Qt[ntq-1]; } } /* * Establish the heap-order property by applying filterDown to * all nodes that have children, ordered from bottom up. */ for (ip = ntq/2; ip > 0; ip--) { filterDown(ip); } return; } /* End of buildQt */ /* *-------------------------------------------------------------------- * * closeR: Close the user-selected polygon R by implicitly adding * an edge from the last vertex to the first. If the added * edge would result in a self-intersecting curve, the * polygon is rejected by setting np = 0. * * Global variables required: * * np = Number of vertices in the boundary of R. * * R = Array dimensioned npmax by 2 containing the Cartesian * coordinates of the sequence of vertices defining the * boundary of R. * * glmesh2 function called: intsec * *-------------------------------------------------------------------- */ static void closeR(void) { GLboolean intr = 0; GLint i; /* * Test for at least three boundary vertices. */ if (np < 3) np = 0; if (np == 0) return; /* * Test the last segment for an intersection. */ for (i = 1; i <= np-3; i++) { if (intsec(R[np-1][0],R[np-1][1],R[0][0],R[0][1], R[i][0],R[i][1],R[i+1][0],R[i+1][1])) { intr = 1; break; } } if (intr) np = 0; return; } /* End of closeR */ /* *-------------------------------------------------------------------- * * coarsenRv: Coarsen R by removing approximately 70 percent of the * vertices, chosen at random. * * glmesh2 functions called: delvtx, initR * *-------------------------------------------------------------------- */ static void coarsenRv(void) { int k, kv, n, nd, nrr, nvr; /* * Store listR if necessary. */ if (newR) initR(); /* * Count the number of vertices nvr and the number of removable * vertices nrr in R, and store their indices in lste. */ nvr = 0; nrr = 0; for (kv = 0; kv < nv; kv++) { if (!listR[kv]) continue; nvr++; if (listn[kv]) continue; lste[nrr] = kv; nrr++; } if (!nrr) { printf("\n coarsenRv: there are no removable vertices in R.\n"); return; } /* * Delete n = max(.7*nvr,nrr) removable vertices. */ n = (int) (0.7*(double) nvr + 0.5); if (n > nrr) n = nrr; nd = 0; while (nd < n && nrr > 0) { k = rand()%nrr; kv = lste[k]; if (delvtx(kv)) { nd++; /* * Vertex nv was relabeled as kv. If index nv is in lste, it must * be changed to kv. */ for (k = 0; k < nrr; k++) { if (lste[k] == nv) lste[k] = kv; } } nrr--; lste[k] = lste[nrr]; } /* * The refinement stack is emptied if a vertex was deleted. */ if (nd) nref = 0; return; } /* End of coarsenRv */ /* *-------------------------------------------------------------------- * * convertRv: Convert boundary vertices in R to removable or non- * removable by deleting their indices from lnrv or * inserting their indices into lnrv, and adjusting the * flags in listn. * * On input: removable is a flag specifying that the boundary * vertices in R are to be converted to removable (TRUE) * or non-removable (FALSE). * * glmesh2 function called: initR * *-------------------------------------------------------------------- */ static void convertRv(unsigned char removable) { int i, kt, kv; /* * Store listR if necessary. */ if (newR) initR(); if (removable) { /* * For each non-removable boundary vertex kv in R, * remove kv from lnrv. */ for (kv = 0; kv < nv; kv++) { if (!listn[kv]) continue; if (!listR[kv]) continue; for (i = 0; i < nn; i++) { if (lnrv[i] == kv) { nn--; lnrv[i] = lnrv[nn]; listn[kv] = 0; break; } } } } else { /* * For each removable boundary vertex kv in R, add kv to lnrv. */ for (kv = 0; kv < nv; kv++) { if (listn[kv]) continue; if (!listR[kv]) continue; kt = lvtx[kv]; if (ltri[kt][0] == kv) i = 2; else if (ltri[kt][1] == kv) i = 0; else i = 1; if (ltri[kt][i+3] >= 0) continue; lnrv[nn] = kv; nn++; listn[kv] = 1; } } return; } /* End of convertRv */ /* *-------------------------------------------------------------------- * * deleteQt: Return the first triangle in priority queue Qt, and * remove it from the queue. * * On input: ptri is a pointer to a struct triangle into which the * first triangle will be copied unless the queue is * empty (ntq = 0). * * Global variables required: ntq, Qtp * * glmesh2 function required: filterDown * * *-------------------------------------------------------------------- */ static void deleteQt(struct triangle *ptri) { int indx; if (!ntq) return; /* * Copy *Qtp[1] into *ptri */ ptri->index_qtp = Qtp[1]->index_qtp; ptri->index_t = Qtp[1]->index_t; ptri->index_v1 = Qtp[1]->index_v1; ptri->index_v2 = Qtp[1]->index_v2; ptri->index_v3 = Qtp[1]->index_v3; ptri->imin = Qtp[1]->imin; ptri->emin = Qtp[1]->emin; /* * Copy Qt[ntq-1] into *Qtp[1], and correct the appropriate pointer * (Qtp value). */ Qtp[1]->index_qtp = Qt[ntq-1].index_qtp; Qtp[1]->index_t = Qt[ntq-1].index_t; Qtp[1]->index_v1 = Qt[ntq-1].index_v1; Qtp[1]->index_v2 = Qt[ntq-1].index_v2; Qtp[1]->index_v3 = Qt[ntq-1].index_v3; Qtp[1]->imin = Qt[ntq-1].imin; Qtp[1]->emin = Qt[ntq-1].emin; indx = Qt[ntq-1].index_qtp; Qtp[indx] = Qtp[1]; /* * Copy Qtp[ntq] into Qtp[1], correct the index_qtp value, and * decrement ntq. */ Qtp[1] = Qtp[ntq]; Qtp[1]->index_qtp = 1; ntq--; /* * Restore the heap-order property by filtering the insertion point * down. */ filterDown(1); return; } /* End of deleteQt */ /* *-------------------------------------------------------------------- * * deltri: Delete a triangle. * * This function removes a triangle from the triangle mesh unless its * removal would delete a non-removable vertex or destroy the * mesh by leaving the union of remaining triangles not strongly * connected. * * On input: kt0 is the index of the triangle to be removed. * * On output: the triangle is removed unless its index is not valid * (in the range 0 to nt-1), or it is not in R, or its * removal would render the union of triangles weakly * connected. The return value is the number of * triangles removed. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static unsigned char deltri(int kt0) { int i, j, ke, kt, kt1, ktn, kv, kvs, nbe; /* * Test for invalid index. */ if (kt0 < 0 || kt0 >= nt) return 0; /* * kt0 cannot be removed if it has a boundary vertex kv which is * not an endpoint of a boundary edge contained in kt0. Also, * test for a vertex kv of kt0 not in R. */ for (i = 0; i < 3; i++) { kv = ltri[kt0][i]; if (!listR[kv]) return 0; kt = lvtx[kv]; if (ltri[kt][0] == kv) j = 2; else if (ltri[kt][1] == kv) j = 0; else j = 1; if (ltri[kt][j+3] < 0 && kt != kt0 && ltri[kt0][(i+1)%3+3] >= 0) return 0; } /* * Test for a non-removable vertex shared by two boundary edges * in kt0. */ for (i = 0; i < 3; i++) { kt = ltri[kt0][(i+1)%3+3]; ktn = ltri[kt0][(i+2)%3+3]; if (kt < 0 && ktn < 0) { kv = ltri[kt0][i]; if (listn[kv]) return 0; } } /* * Save a vertex index for updating lcrv if necessary. */ kvs = kv; /* * Loop on vertices kv and opposite edges ke of triangle kt0, * accumulating nbe (number of boundary edges in kt0). */ nbe = 0; for (i = 0; i < 3; i++) { kv = ltri[kt0][i]; /* * Alter lvtx[kv] if necessary, to reflect the new boundary curve. */ kt = ltri[kt0][(i+1)%3+3]; if (kt >= 0) { lvtx[kv] = kt; } else { if (ltri[kt0][(i+2)%3+3] < 0) { /* * Remove vertex kv, and relabel vertex nv-1 to kv. */ nv--; vtxy[kv][0] = vtxy[nv][0]; vtxy[kv][1] = vtxy[nv][1]; uv[kv][0] = uv[nv][0]; uv[kv][1] = uv[nv][1]; kt = lvtx[nv]; lvtx[kv] = kt; kt1 = kt; do { if (ltri[kt][0] == nv) j = 0; else if (ltri[kt][1] == nv) j = 1; else j = 2; ltri[kt][j] = kv; j = (j+1)%3; kt = ltri[kt][j+3]; } while (kt >= 0 && kt != kt1); /* * Correct lcrv if necessary. */ for (j = 0; j < nc; j++) { if (lcrv[j] == kv) lcrv[j] = ltri[kt0][(i+1)%3]; if (lcrv[j] == nv) lcrv[j] = kv; } /* * Correct lnrv and listn if necessary. */ for (j = 0; j < nn; j++) { if (lnrv[j] == nv) { lnrv[j] = kv; listn[kv] = 1; break; } } /* * Adjust listR. */ listR[kv] = listR[nv]; } } ktn = ltri[kt0][i+3]; ke = ltri[kt0][i+6]; if (ktn >= 0) { /* * Remove kt0 as a neighbor of ktn, and alter ledg[ke] if necessary. */ if (ltri[ktn][3] == kt0) j = 0; else if (ltri[ktn][4] == kt0) j = 1; else j = 2; ltri[ktn][j+3] = -1; ledg[ke] = ktn; } else { /* * Remove boundary edge ke, and relabel edge ne-1 to ke. */ nbe++; ne--; if (ke != ne) { kt = ledg[ne]; ledg[ke] = kt; if (ltri[kt][6] == ne) j = 6; else if (ltri[kt][7] == ne) j = 7; else j = 8; ltri[kt][j] = ke; kt = ltri[kt][j-3]; if (kt >= 0) { if (ltri[kt][6] == ne) j = 6; else if (ltri[kt][7] == ne) j = 7; else j = 8; ltri[kt][j] = ke; } } } } /* * Remove triangle kt0, and relabel triangle nt-1 to kt0. */ nt--; if (kt0 != nt) { for (j = 0; j < 9; j++) { ltri[kt0][j] = ltri[nt][j]; } for (i = 0; i < 3; i++) { kv = ltri[nt][i]; if (lvtx[kv] == nt) lvtx[kv] = kt0; kt = ltri[nt][i+3]; if (kt >= 0) { if (ltri[kt][3] == nt) ltri[kt][3] = kt0; else if (ltri[kt][4] == nt) ltri[kt][4] = kt0; else if (ltri[kt][5] == nt) ltri[kt][5] = kt0; } ke = ltri[nt][i+6]; if (ledg[ke] == nt) { ledg[ke] = (kt > kt0) ? kt : kt0; } } } /* * Adjust counts nb and nc, and update lcrv if a new curve was added. */ if (nbe == 0) { lcrv[nc] = kvs; nc++; nb += 3; } else if (nbe == 1) { nb++; } else if (nbe == 2) { nb--; } else { nc--; nb -= 3; } return 1; } /* End of deltri */ /* *-------------------------------------------------------------------- * * delvtx: Delete a vertex. * * All possible edge swaps (Function swap) are applied to the edges * incident on vertex kv0 so that kv0, if an interior vertex, has * three neighbors, and then kv0 and the edges containing kv0 are * removed, replacing three triangles by one triangle in the case of * an interior vertex. The vertex is thus removed without leaving a * hole. * * This function, called with kv0 = nv-1, reverses the effect of a * call to Function insertv. * * On input: kv0 is the index of the vertex to be removed. * * On output: The return value is the number of vertices removed * (0 or 1). If index kv0 is invalid, or kv0 is an * element of lnrv, no variables are altered. It is * possible that kv0 indexes an interior vertex but the * number of neighbors of kv0 could not be reduced to * three due to collinear vertices and floating point * inaccuracy. Also, it is possible that a boundary * vertex cannot be removed without destroying the * triangle mesh. The mesh may have been altered by * edge swaps but kv0 is not removed in these cases. * * Note that vertex kv0 may be deleted even if it is not contained * in R. * * glmesh2 function called: swap * *-------------------------------------------------------------------- */ static unsigned char delvtx(int kv0) { unsigned char bndry; int i, i0, iter, j, ke, ke0, ken, kt, kt0, kt1, ktn, kv, kv1, kv2, kv3, kvs = 0, nnb, swp; /* * Test for invalid input. */ if (kv0 < 0 || kv0 >= nv || nv < 4) return 0; /* * Test for kv0 non-removable. */ if (listn[kv0]) return 0; /* * Count the number of neighbors of kv0. * * kt1 = index of first triangle containing kv0. * nnb = number of neighbors of kv0. */ kt1 = lvtx[kv0]; kt = kt1; nnb = 0; do { if (ltri[kt][0] == kv0) i = 1; else if (ltri[kt][1] == kv0) i = 2; else i = 0; kt = ltri[kt][i+3]; nnb++; } while (kt >= 0 && kt != kt1); bndry = (kt < 0); if (bndry) nnb++; /* * Loop on edges ke incident on kv0 and shared by triangles * kt and ktn, swapping out ke if possible. * * bndry = 1 iff kv0 is a boundary vertex. * kt1 = index of first triangle containing kv0. * nnb = number of neighbors of kv0. * swp = flag with value 1 initially and following a swap. */ kt = kt1; swp = 1; while ((swp || kt != kt1) && (bndry || nnb > 3)) { if (ltri[kt][0] == kv0) i = 1; else if (ltri[kt][1] == kv0) i = 2; else i = 0; ktn = ltri[kt][i+3]; if (ktn < 0) break; ke = ltri[kt][i+6]; kv2 = ltri[kt][i]; kv1 = ltri[kt][(i+1)%3]; if (ltri[ktn][0] == kv0) i = 2; else if (ltri[ktn][1] == kv0) i = 0; else i = 1; kv3 = ltri[ktn][i]; /* * Edge ke can be swapped iff kv1 strictly Left kv3 -> kv2 and * kv0 Left kv2 -> kv3. */ swp = ((vtxy[kv2][0]-vtxy[kv3][0])*(vtxy[kv1][1]-vtxy[kv3][1]) > (vtxy[kv2][1]-vtxy[kv3][1])*(vtxy[kv1][0]-vtxy[kv3][0]) && (vtxy[kv3][0]-vtxy[kv2][0])*(vtxy[kv0][1]-vtxy[kv2][1]) >= (vtxy[kv3][1]-vtxy[kv2][1])*(vtxy[kv0][0]-vtxy[kv2][0])); if (swp) { if (swap(ke)) { nnb--; if (kt == kt1) kt1 = ktn; } else { swp = 0; } } kt = ktn; } /* * Test for failure with kv0 interior and having more than three * neighbors, or kv0 a boundary vertex with a boundary vertex as * a neighbor other than the first or last. */ if (!bndry) { if (nnb > 3) return 0; } else { kt = lvtx[kv0]; while (kt >= 0) { if (ltri[kt][0] == kv0) i = 1; else if (ltri[kt][1] == kv0) i = 2; else i = 0; ktn = ltri[kt][i+3]; if (ktn < 0) break; i = (i+1)%3; kv1 = ltri[kt][i]; kt1 = lvtx[kv1]; if (ltri[kt1][0] == kv1) j = 2; else if (ltri[kt1][1] == kv1) j = 0; else j = 1; if (ltri[kt1][j+3] < 0) return 0; kt = ktn; } } /* * Set kt0 and ke0 to the indices of the first triangle and edge, * respectively, to be deleted, and set i0 to the ltri column * index of kv0 in kt1. */ kt1 = lvtx[kv0]; if (ltri[kt1][0] == kv0) i0 = 0; else if (ltri[kt1][1] == kv0) i0 = 1; else i0 = 2; i = i0; if (!bndry) { i = (i+1)%3; ke0 = ltri[kt1][i+6]; kt0 = ltri[kt1][i+3]; j = (i+1)%3; kv1 = ltri[kt1][j]; if (lvtx[kv1] == kt0) lvtx[kv1] = kt1; } else { /* * Save a neighbor kvs of kv0 for updating lcrv if necessary. */ kvs = ltri[kt1][(i0+1)%3]; i = (i+2)%3; ke0 = ltri[kt1][i+6]; kt0 = kt1; kt1 = -1; } /* * Remove edge ke0, and relabel edge ne-1 to ke0. */ ne--; if (ke0 != ne) { kt = ledg[ne]; ledg[ke0] = kt; if (ltri[kt][6] == ne) j = 6; else if (ltri[kt][7] == ne) j = 7; else j = 8; ltri[kt][j] = ke0; kt = ltri[kt][j-3]; if (kt >= 0) { if (ltri[kt][6] == ne) j = 6; else if (ltri[kt][7] == ne) j = 7; else j = 8; ltri[kt][j] = ke0; } } /* * Loop on triangles kt0 with triangle ktn and edge ken opposite * vertex kv0. Note that kt1 = -1 iff bndry = 1. */ for (iter = 1; iter < nnb; iter++) { if (ltri[kt0][0] == kv0) i = 0; else if (ltri[kt0][1] == kv0) i = 1; else i = 2; ktn = ltri[kt0][i+3]; ken = ltri[kt0][i+6]; if (ktn >= 0) { if (ltri[ktn][3] == kt0) ltri[ktn][3] = kt1; else if (ltri[ktn][4] == kt0) ltri[ktn][4] = kt1; else ltri[ktn][5] = kt1; } if (!bndry) { j = (i+1)%3; kv1 = ltri[kt0][j]; j = (i+2)%3; kv2 = ltri[kt0][j]; if (lvtx[kv1] == kt0) lvtx[kv1] = kt1; if (lvtx[kv2] == kt0) lvtx[kv2] = kt1; ledg[ken] = (kt1 > ktn) ? kt1 : ktn; ltri[kt1][(i0+iter)%3+3] = ktn; ltri[kt1][(i0+iter)%3+6] = ken; } else { j = (i+2)%3; kv2 = ltri[kt0][j]; lvtx[kv2] = ktn; ledg[ken] = ktn; } /* * Set ke0 to the edge of triangle kt0 that contains kv0 and was not * removed at the previous iteration (or before the loop), and set * ktn to the next triangle to be processed. */ i = (i+1)%3; ke0 = ltri[kt0][i+6]; ktn = ltri[kt0][i+3]; /* * Remove edge ke0, and relabel edge ne-1 to ke0. */ ne--; if (ke0 != ne) { kt = ledg[ne]; ledg[ke0] = kt; if (ltri[kt][6] == ne) j = 6; else if (ltri[kt][7] == ne) j = 7; else j = 8; ltri[kt][j] = ke0; kt = ltri[kt][j-3]; if (kt >= 0) { if (ltri[kt][6] == ne) j = 6; else if (ltri[kt][7] == ne) j = 7; else j = 8; ltri[kt][j] = ke0; } } /* * Store the new first vertex of kt1 if kv0 is interior. */ if (!bndry && iter == 1) ltri[kt1][i0] = ltri[kt0][(i+1)%3]; /* * Remove triangle kt0, and relabel triangle nt-1 to kt0. */ nt--; if (kt0 != nt) { if (kt1 == nt) kt1 = kt0; if (ktn == nt) ktn = kt0; for (j = 0; j < 9; j++) { ltri[kt0][j] = ltri[nt][j]; } for (i = 0; i < 3; i++) { kv = ltri[nt][i]; if (lvtx[kv] == nt) lvtx[kv] = kt0; kt = ltri[nt][i+3]; if (kt >= 0) { if (ltri[kt][3] == nt) ltri[kt][3] = kt0; else if (ltri[kt][4] == nt) ltri[kt][4] = kt0; else if (ltri[kt][5] == nt) ltri[kt][5] = kt0; } ke = ltri[nt][i+6]; if (ledg[ke] == nt) { ledg[ke] = (kt > kt0) ? kt : kt0; } } } /* * Bottom of for loop. */ kt0 = ktn; } /* * Remove vertex kv0 from vtxy, uv, and lvtx, and relabel vertex * nv-1 to kv0. */ nv--; if (bndry) nb += nnb-3; for (j = 0; j < nc; j++) { if (lcrv[j] == kv0) lcrv[j] = kvs; if (lcrv[j] == nv) lcrv[j] = kv0; } if (kv0 != nv) { vtxy[kv0][0] = vtxy[nv][0]; vtxy[kv0][1] = vtxy[nv][1]; uv[kv0][0] = uv[nv][0]; uv[kv0][1] = uv[nv][1]; kt = lvtx[nv]; lvtx[kv0] = kt; kt1 = kt; do { if (ltri[kt][0] == nv) j = 0; else if (ltri[kt][1] == nv) j = 1; else j = 2; ltri[kt][j] = kv0; j = (j+1)%3; kt = ltri[kt][j+3]; } while (kt >= 0 && kt != kt1); if (listn[nv]) { /* * Relabel nv in lnrv. */ for (i = 0; i < nn; i++) { if (lnrv[i] == nv) { lnrv[i] = kv0; break; } } } } /* * Adjust listn and listR. */ listn[kv0] = listn[nv]; listR[kv0] = listR[nv]; return 1; } /* End of delvtx */ /* *-------------------------------------------------------------------- * * dlgDisplay: Dialog window display callback. * * Global variables required: * * dlg_r0 = Rectangle representing 'No' button. * * dlg_r1 = Rectangle representing 'Yes' button. * * dlg_string = Text string displayed in dialog box window. * * glmesh2 functions called: displayString, drawRect * *-------------------------------------------------------------------- */ static void dlgDisplay(void) { GLfloat black[3] = {0.0, 0.0, 0.0}; /* Border colors */ GLfloat white[3] = {1.0, 1.0, 1.0}; int fontno = 3; /* Font number for strings (24-point proportional spaced Times-Roman) */ GLfloat pos[3]; /* Text string position */ char text[7]; /* Text string */ int x0, y0; /* Coordinates of the lower left corner of a rectangle */ /* * Clear the color buffer to the background color. */ glClear(GL_COLOR_BUFFER_BIT); /* * Display dlg_string. */ pos[0] = 75.0f; pos[1] = 160.0f; pos[2] = 0.0f; displayString(pos, dlg_string, fontno, white); /* * Draw buttons at the bottom of the window. */ x0 = dlg_r1.x0; y0 = dlg_r1.y0; drawRect(x0, y0, dlg_r1.w, dlg_r1.h, black, white); pos[0] = (GLfloat) (x0+7); pos[1] = (GLfloat) (y0+7); sprintf(text, " YES "); displayString(pos, text, fontno, white); x0 = dlg_r0.x0; y0 = dlg_r0.y0; drawRect(x0, y0, dlg_r0.w, dlg_r0.h, black, white); pos[0] = (GLfloat) (x0+7); pos[1] = (GLfloat) (y0+7); sprintf(text, " NO "); displayString(pos, text, fontno, white); glFlush(); return; } /* End of dlgDisplay */ /* *-------------------------------------------------------------------- * * dlgKey: Keyboard callback associated with the dialog window: * Enter ==> Yes * Escape ==> No * * The mouse coordinates (x,y) are not used. * * Global variables required: * * dlg_r1 = Rectangle representing 'Yes' button. * * dlg_yes = Flag with value defined by the 'Yes' button: set * to FALSE by the Escape key. * * glmesh2 function called: dlgMouse * *-------------------------------------------------------------------- */ static void dlgKey(unsigned char key, int x, int y) { if (key == 10 || key == 13) { /* Enter key pressed */ dlgMouse(GLUT_LEFT_BUTTON, GLUT_UP, dlg_r1.x0+1, 198-dlg_r1.y0); } else if (key == 27) { /* Escape key pressed */ dlg_yes = GL_FALSE; glutHideWindow(); glutSetWindow(window[0]); glutShowWindow(); /* Restore visibility of window 0 */ } return; } /* End of dlgKey */ /* *-------------------------------------------------------------------- * * dlgMouse: Callback triggered by a mouse button press or release * with the mouse position in the dialog window. * * If the left mouse button is pressed with the mouse position * inside the 'Yes' or 'No' rectangle, the appearance of the * rectangle is altered. If the left button is released inside one * of these rectangles, dlg_yes is set to TRUE or FALSE, and the * window state is changed to hidden. * * Global variables required: * * dlg_exit = Flag with value TRUE iff the program is to be * terminated when the user selects the 'Yes' button. * * dlg_r0 = Rectangle representing 'No' button. * * dlg_r1 = Rectangle representing 'Yes' button. * * dlg_yes = Flag with value TRUE iff the user selects the 'Yes' * button in the dialog box. * * fname = File name for output. * * Triangulation variables: ltri, nb, nc, ne, nt, nv, vtxy, uv * * glmesh2 functions called: drawRect, trwrite * *-------------------------------------------------------------------- */ static void dlgMouse(int button, int state, int x, int y) { GLfloat black[3] = {0.0, 0.0, 0.0}; /* Border colors */ GLfloat white[3] = {1.0, 1.0, 1.0}; int x0, y0; /* Lower left corner of a rectangle */ int x1, y1; /* Upper right corner of a rectangle */ y = 199 - y; /* Map glut y-coordinate to viewport */ /* * Left button press or release. */ x0 = dlg_r1.x0; y0 = dlg_r1.y0; x1 = x0 + dlg_r1.w - 1; y1 = y0 + dlg_r1.h - 1; if (x > x0 && x < x1 && y > y0 && y < y1) { /* * 'Yes' button. */ if (state == GLUT_DOWN) { drawRect(x0, y0, dlg_r1.w, dlg_r1.h, white, black); glFlush(); } else { /* 'Yes' button released */ if (dlg_exit) { /* * Write the triangle mesh to a file, and terminate execution. */ trwrite(fname, nc, nb, nv, nt, ne, vtxy, uv, ltri); exit(0); } dlg_yes = GL_TRUE; glutHideWindow(); glutSetWindow(window[0]); glutShowWindow(); } return; } x0 = dlg_r0.x0; y0 = dlg_r0.y0; x1 = x0 + dlg_r0.w - 1; y1 = y0 + dlg_r0.h - 1; if (x > x0 && x < x1 && y > y0 && y < y1) { /* * 'No' button. */ if (state == GLUT_DOWN) { drawRect(x0, y0, dlg_r0.w, dlg_r0.h, white, black); glFlush(); } else { /* 'No' button released */ dlg_yes = GL_FALSE; glutHideWindow(); glutSetWindow(window[0]); glutShowWindow(); } return; } return; } /* End of dlgMouse */ /* *-------------------------------------------------------------------- * * dlgReshape: Set the viewport, projection matrix, and modelview * matrix for the dialog window (window 1): * two-dimensional rendering with a 400 by 200 viewport. * * This function is called before the first call to dlgDisplay. * * If the window is too small, a request for a larger window is * executed. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void dlgReshape(int w, int h) { if (w < 400 || h < 200) { /* Test for sufficient */ glutReshapeWindow(400, 200); /* window extent. */ return; } /* * Set the viewport to the entire window (the default). */ glViewport(0, 0, (GLint) w, (GLint) h); /* * Set up an orthogonal projection for data space equal to * viewport (identity). */ glMatrixMode(GL_PROJECTION); glLoadIdentity(); gluOrtho2D(0.0, (GLdouble) w, 0.0, (GLdouble) h); /* * Set the modelview matrix to the identity. */ glMatrixMode(GL_MODELVIEW); glLoadIdentity(); return; } /* End of dlgReshape */ /* *-------------------------------------------------------------------- * * display: Clear the color buffer and draw the triangle mesh edges * with optional color-fill and index annotation of ver- * tices, edges, and triangles. The title, and selected * polygon R are also drawn, and a bounding box may be * rendered. * * glmesh2 functions called: displayString, drawCircle, trqual1 * *-------------------------------------------------------------------- */ static void display(void) { unsigned int fontno; /* Font for labels */ GLfloat pos[3]; /* Label position */ char text[30]; /* Label text */ GLfloat red[3] = {1.0, 0.0, 0.0}; GLfloat green[3] = {0.0, 1.0, 0.0}; GLfloat blue[3] = {0.0, 0.0, 1.0}; GLfloat grey[3] = {0.7, 0.7, 0.7}; double q; int i; GLdouble ch; GLint kc, kc1, kc2, kc3, ke, kt, kt1, kt2, kt3, kv, kv1, kv2, kv3; /* * Clear the color buffer, and set line and point widths. */ glClear(GL_COLOR_BUFFER_BIT); glLineWidth(line_width); glPointSize(line_width+1); if (color_fill) glPolygonMode(GL_FRONT_AND_BACK, GL_FILL); if (color_fill == 1 || color_fill == 2) { /* * Color-fill the triangles (only those in R if color_fill = 2) * using four colors (from array colors) such that adjacent * triangles are rendered in different colors. Color indices, * along with flag value -1 (index not yet assigned) are stored * in lste. */ for (kt = 0; kt < nt; kt++) lste[kt] = -1; for (kt = 0; kt < nt; kt++) { kv1 = ltri[kt][0]; kv2 = ltri[kt][1]; kv3 = ltri[kt][2]; if (color_fill == 2) { if (!listR[kv1] || !listR[kv2] || !listR[kv3]) continue; } kt1 = ltri[kt][3]; kt2 = ltri[kt][4]; kt3 = ltri[kt][5]; kc1 = -1; kc2 = -1; kc3 = -1; if (kt1 >= 0) kc1 = lste[kt1]; if (kt2 >= 0) kc2 = lste[kt2]; if (kt3 >= 0) kc3 = lste[kt3]; kc = 0; while (kc == kc1 || kc == kc2 || kc == kc3) { kc = (kc+1)%4; } lste[kt] = kc; glColor3fv((GLfloat *) &colors[kc]); glBegin(GL_TRIANGLES); glVertex2d(vtxy[kv1][0], vtxy[kv1][1]); glVertex2d(vtxy[kv2][0], vtxy[kv2][1]); glVertex2d(vtxy[kv3][0], vtxy[kv3][1]); glEnd(); } } if (color_fill == 3) { /* * Render constant-color triangles using color-coding based on * quality. */ for (kt = 0; kt < nt; kt++) { kv1 = ltri[kt][0]; kv2 = ltri[kt][1]; kv3 = ltri[kt][2]; q = trqual1(vtxy[kv1][0], vtxy[kv1][1], vtxy[kv2][0], vtxy[kv2][1], vtxy[kv3][0], vtxy[kv3][1]); if (q < 0.600) glColor3fv(red); else if (q < 0.866) glColor3fv(blue); else glColor3fv(green); glBegin(GL_TRIANGLES); glVertex2d(vtxy[kv1][0], vtxy[kv1][1]); glVertex2d(vtxy[kv2][0], vtxy[kv2][1]); glVertex2d(vtxy[kv3][0], vtxy[kv3][1]); glEnd(); } } if (color_fill != 1) { /* * Draw the triangle mesh edges. */ glPolygonMode(GL_FRONT_AND_BACK, GL_LINE); glColor3fv(color_e); for (ke = 0; ke < ne; ke++) { kt = ledg[ke]; if (ltri[kt][6] == ke) i = 1; else if (ltri[kt][7] == ke) i = 2; else i = 0; kv1 = ltri[kt][i]; kv2 = ltri[kt][(i+1)%3]; glBegin(GL_LINES); glVertex2d(vtxy[kv1][0], vtxy[kv1][1]); glVertex2d(vtxy[kv2][0], vtxy[kv2][1]); glEnd(); } } /* * Display the selected polygon R (both line segments and vertices). */ glColor3fv(red); for (i = 0; i < np; i++) { glBegin(GL_LINES); glVertex2d(R[i][0], R[i][1]); glVertex2d(R[(i+1)%np][0], R[(i+1)%np][1]); glEnd(); glBegin(GL_POINTS); glVertex2d(R[i][0], R[i][1]); glEnd(); } if (cr > 0.0) { /* * Draw a circle of radius cr centered at (cx,cy). */ glColor3fv(grey); drawCircle(); } if (plot_v) { /* * Display vertex indices. */ fontno = ind_font; pos[2] = 0.0; for (kv = 0; kv < nv; kv++) { pos[0] = (GLfloat) (vtxy[kv][0] + .005*dx); pos[1] = (GLfloat) (vtxy[kv][1] + .005*dy); sprintf(text, "%u", kv); displayString(pos, text, fontno, color_t); } } if (plot_nrv) { /* * Display non-removable vertex indices. */ fontno = ind_font; pos[2] = 0.0; for (i = 0; i < nn; i++) { kv = lnrv[i]; pos[0] = (GLfloat) (vtxy[kv][0] + .005*dx); pos[1] = (GLfloat) (vtxy[kv][1] + .005*dy); sprintf(text, "%u", kv); displayString(pos, text, fontno, green); } } if (plot_e) { /* * Display edge indices. */ fontno = ind_font; pos[2] = 0.0; for (ke = 0; ke < ne; ke++) { kt = ledg[ke]; if (ltri[kt][6] == ke) i = 1; else if (ltri[kt][7] == ke) i = 2; else i = 0; kv1 = ltri[kt][i]; kv2 = ltri[kt][(i+1)%3]; pos[0] = (GLfloat) ((vtxy[kv1][0]+vtxy[kv2][0])/2.0 + .005*dx); pos[1] = (GLfloat) ((vtxy[kv1][1]+vtxy[kv2][1])/2.0 + .005*dy); sprintf(text, "%u", ke); displayString(pos, text, fontno, color_t); } } if (plot_t) { /* * Display triangle indices. */ fontno = ind_font; pos[2] = 0.0; for (kt = 0; kt < nt; kt++) { kv1 = ltri[kt][0]; kv2 = ltri[kt][1]; kv3 = ltri[kt][2]; pos[0] = (GLfloat) ((vtxy[kv1][0]+vtxy[kv2][0]+vtxy[kv3][0]) /3.0); pos[1] = (GLfloat) ((vtxy[kv1][1]+vtxy[kv2][1]+vtxy[kv3][1]) /3.0); sprintf(text, "%u", kt); displayString(pos, text, fontno, color_t); } } if (plot_box) { /* * Display the bounding box with the lower left and upper right * corners labeled. */ glColor3fv(grey); fontno = 5; glBegin(GL_LINE_LOOP); glVertex2d(xmin, ymin); glVertex2d(xmax, ymin); glVertex2d(xmax, ymax); glVertex2d(xmin, ymax); glEnd(); /* * The character height ch is approximately 11 pixels (converted * to object coordinates). */ ch = (sx <= sy) ? sx : sy; /* ch = min(sx, sy) */ ch = 11.0*ch; pos[0] = (GLfloat) (xmin); pos[1] = (GLfloat) (ymin - 1.5*ch); pos[2] = 0.0; sprintf(text, "(%3.2f, %3.2f)", xmin, ymin); displayString(pos, text, fontno, color_t); pos[0] = (GLfloat) (xmax - 6.0*ch); pos[1] = (GLfloat) (ymax + .5*ch); sprintf(text, "(%3.2f, %3.2f)", xmax, ymax); displayString(pos, text, fontno, color_t); } if (*title != '\0') { /* * Display the title string */ displayString(title_pos, title, title_font, color_t); } glColor3fv(color_e); glutSwapBuffers(); return; } /* End of display */ /* *-------------------------------------------------------------------- * * displayPosition: Display the coordinates of a point, such as the * current mouse position, in the upper left corner * of the viewport. * * On input: x and y are the coordinates to be displayed. * * Global variables required: color_b, color_t, dx, dy, sx, sy, * xc, yc * * glmesh2 functions called: displayString * *-------------------------------------------------------------------- */ static void displayPosition(GLdouble x, GLdouble y) { static GLdouble xsav = 0.0, ysav = 0.0; /* Previous (x,y) */ unsigned int fontno = 5; /* Font for strings */ GLfloat pos[3]; /* Label position */ char text[30]; /* Label text */ GLdouble ch; /* Character height */ GLboolean xor; /* * The character height ch is approximately 11 pixels (converted * to object coordinates). The strings written on the previous * call are erased by writing them in the background color. */ ch = (sx <= sy) ? sx : sy; /* ch = min(sx, sy) */ ch = 11.0*ch; pos[0] = (GLfloat) (xc-0.5*dx+ch); pos[1] = (GLfloat) (yc+0.5*dy-1.5*ch); pos[2] = 0.0; /* * Strings are written in overwrite mode. The strings written on * the previous call are erased by writing them in the background * color. */ xor = glIsEnabled(GL_COLOR_LOGIC_OP); glDisable(GL_COLOR_LOGIC_OP); sprintf(text, "x = %.3e", xsav); displayString(pos, text, fontno, color_b); sprintf(text, "x = %.3e", x); displayString(pos, text, fontno, color_t); pos[1] = (GLfloat) (yc+0.5*dy-3.0*ch); sprintf(text, "y = %.3e", ysav); displayString(pos, text, fontno, color_b); sprintf(text, "y = %.3e", y); displayString(pos, text, fontno, color_t); /* * Update saved (x,y), and restore the write mode. */ xsav = x; ysav = y; if (xor) { glEnable(GL_COLOR_LOGIC_OP); glLogicOp(GL_XOR); } glutSwapBuffers(); return; } /* End of displayPosition */ /* *-------------------------------------------------------------------- * * displayString: Write string to the screen with lower left corner * at raster position pos using font number fontno * (defined in the table below) and color scolor. * * Note that the raster position is actually a location in the object * coordinate space; i.e., the current modelview and projection * matrices are applied to pos in order to obtain the screen coordi- * nates of the first character. If those coordinates are not in the * viewing volume, nothing is written. Otherwise the string is * written but clipped against the window boundaries. * * If fontno is invalid, an error message is written to the standard * output unit and the return value is 1. Otherwise the return value * is 0. * * Global variables required: dx, dy * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static int displayString(GLfloat pos[3], char *string, unsigned int fontno, GLfloat scolor[3]) { unsigned int i; /* string index */ unsigned int nfonts = 9; /* Number of fonts */ GLfloat sf; /* Scale factor for stroked fonts */ static void *font[9] = { GLUT_BITMAP_8_BY_13, /* 8 x 13 fixed width */ GLUT_BITMAP_9_BY_15, /* 9 x 15 fixed width */ GLUT_BITMAP_TIMES_ROMAN_10, /* 10-point proportional spaced */ GLUT_BITMAP_TIMES_ROMAN_24, /* 24-point proportional spaced */ GLUT_BITMAP_HELVETICA_10, /* 10-point proportional spaced */ GLUT_BITMAP_HELVETICA_12, /* 12-point proportional spaced */ GLUT_BITMAP_HELVETICA_18, /* 18-point proportional spaced */ GLUT_STROKE_ROMAN, /* proportional spaced stroke */ GLUT_STROKE_MONO_ROMAN /* mono-spaced Roman simplex */ }; if (fontno >= nfonts) { printf("displayString: Invalid font number = %u\n", fontno); return 1; } glColor3fv(scolor); if (fontno < 7 ) { glRasterPos3fv(pos); i = 0; while (string[i] != '\0') { glutBitmapCharacter(font[fontno], string[i]); i = i + 1; } } else { sf = 0.0005*(dx+dy)/4.0; glLineWidth(3.0); glPushMatrix(); glTranslatef(pos[0], pos[1], pos[2]); glScalef(sf, sf, sf); i = 0; while (string[i] != '\0') { glutStrokeCharacter(font[fontno], string[i]); i = i + 1; } glPopMatrix(); glLineWidth(line_width); } return 0; } /* End of displayString */ /* *-------------------------------------------------------------------- * * dragCorner: Mouse callback (callback triggered by a mouse button * press or release with the mouse position in the * primary window) used to initiate an alteration of the * position of a non-removable vertex. * * A left button press may be used to select and drag a non-removable * vertex to any location from which its neighbors are visible. * * glmesh2 functions called: bbox, intsec, xlate * *-------------------------------------------------------------------- */ static void dragCorner(int button, int state, int x, int y) { GLdouble wx, wy, wz; /* Computed object coordinates */ GLboolean intr, visible; GLdouble tols, x1, x2, x3, y1, y2, y3; GLint i, kt, kv, kv1, kv2, kv3 = 0, kvn, nnb; if (button != GLUT_LEFT_BUTTON) return; if (state == GLUT_DOWN) { /* * Compute the object coordinates (wx,wy) of the mouse cursor. */ xlate(x,y,0.0, &wx,&wy,&wz); /* * Find a non-removable vertex kv1, if any, within squared distance * tols of the mouse position (wx,wy). */ tols = (sx <= sy) ? sx : sy; /* tols = min(sx, sy) */ tols = 25.0*tols*tols; /* 5-pixel tolerance */ intr = 0; for (i = 0; i < nn; i++) { kv1 = lnrv[i]; if ((intr = ((wx-vtxy[kv1][0])*(wx-vtxy[kv1][0]) + (wy-vtxy[kv1][1])*(wy-vtxy[kv1][1]) <= tols))) break; } if (!intr) return; /* * Set up for dragging in Function motion: store the coordinates of * kv1 in (xv1,yv1), store the indices of the neighboring vertices * of kv1 in the first nnb positions of lste, followed by kv1 in * position nnb, set XOR write mode, set mode = ALTER_DOMAIN, and * set nnbv = nnb. */ xv1 = vtxy[kv1][0]; yv1 = vtxy[kv1][1]; nnb = 0; kt = lvtx[kv1]; do { if (ltri[kt][0] == kv1) i = 1; else if (ltri[kt][1] == kv1) i = 2; else i = 0; lste[nnb] = ltri[kt][i]; nnb++; kvn = ltri[kt][(i+1)%3]; kt = ltri[kt][i+3]; } while (kt >= 0); lste[nnb] = kvn; nnb++; lste[nnb] = kv1; mode = ALTER_DOMAIN; nnbv = nnb; glEnable(GL_COLOR_LOGIC_OP); glLogicOp(GL_XOR); } else { /* * Left button release. */ if (nnbv < 2) return; nnb = nnbv; kv1 = lste[nnb]; /* * The new point p = (xv1,yv1) is visible iff p Left kv2 -> kv3 * for all pairs of adjacent neighbors (kv2,kv3), and there is no * intersection of the path p-kv1 with a boundary edge not containing * kv1 in the same boundary curve. */ visible = 1; for (i = 0; i < nnb-1; i++) { kv2 = lste[i]; kv3 = lste[i+1]; x2 = vtxy[kv2][0]; y2 = vtxy[kv2][1]; x3 = vtxy[kv3][0]; y3 = vtxy[kv3][1]; if ( (x3-x2)*(yv1-y2) < (y3-y2)*(xv1-x2) ) { visible = 0; break; } } if (visible) { /* * Set kv2 and kv3 to the first and last neighbors of kv1, and * test for intersections. */ intr = 0; kv2 = lste[0]; x1 = vtxy[kv1][0]; y1 = vtxy[kv1][1]; kv = kv2; while (kv != kv3) { kt = lvtx[kv]; if (ltri[kt][0] == kv) i = 1; else if (ltri[kt][1] == kv) i = 2; else i = 0; kvn = ltri[kt][i]; if (intsec(xv1,yv1,x1,y1,vtxy[kv][0],vtxy[kv][1], vtxy[kvn][0],vtxy[kvn][1])) { intr = 1; break; } kv = kvn; } if (intr) visible = 0; } if (visible) { /* * Move kv1 to its new position, and update the bounding box. */ vtxy[kv1][0] = xv1; vtxy[kv1][1] = yv1; bbox(); } /* * Set nnbv = 0 indicating that no vertex is currently being dragged. */ nnbv = 0; glutPostRedisplay(); } return; } /* End of dragCorner */ /* *-------------------------------------------------------------------- * * drawCircle: Draw a circle of radius cr centered at (cx,cy), * unless cr is too small (less than 5 pixels). * * Global variables required: cr, cx, cy, sx, sy * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void drawCircle(void) { GLdouble c, da, eps, r, s, x, xs, y; /* * Compute rotation parameters c = cos(da) and s = sin(da) for * angle da = 5/r, where r is the radius in pixels. */ s = (sx >= sy) ? sx : sy; /* da = max(sx, sy) */ r = cr/s; if (r < 5.0) return; da = 5.0/r; c = cos(da); s = sin(da); eps = cr*(1.0-c)/2.0; /* * Loop on points uniformly distributed by angle. */ glBegin(GL_LINE_LOOP); x = cx + cr; y = cy; glVertex2d(x, y); do { xs = x; x = cx + c*(x-cx) - s*(y-cy); y = cy + s*(xs-cx) + c*(y-cy); glVertex2d(x, y); } while (x-cx < cr-eps); glEnd(); glutSwapBuffers(); return; } /* End of drawCircle */ /* *-------------------------------------------------------------------- * * drawRect: Draw the outline of an axis-aligned rectangle in the * z = 0 plane with lower left corner at (x0,y0), width w, * and height h. The origin of the coordinate system is * at the lower left. * * The bottom (lower) and right edges are drawn with color brcolor; * the top and left edges are drawn with color tlcolor, where brcolor * and tlcolor are arrays of length three containing red, green, and * blue components normalized to [0,1]. Taking brcolor and tlcolor * to be black (0,0,0) and white (1,1,1), respectively, gives the * appearance of a button sticking out of the screen. Reversing * those colors makes the button appear to extend inward as if it * were pressed. * * The outline is drawn with lines two pixels wide. If the rectangle * is to be filled, it should be done before the outline is drawn (or * only the interior should be filled) in order to preserve the * outline colors. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void drawRect(GLuint x0, GLuint y0, GLuint w, GLuint h, GLfloat brcolor[3], GLfloat tlcolor[3]) { glLineWidth(2.0); glBegin(GL_LINES); glColor3fv(brcolor); glVertex2i(x0, y0); /* Bottom */ glVertex2i(x0+w-1, y0); glVertex2i(x0+w-1, y0); /* Right */ glVertex2i(x0+w-1, y0+h-1); glColor3fv(tlcolor); glVertex2i(x0+w-1, y0+h-1); /* Top */ glVertex2i(x0, y0+h-1); glVertex2i(x0, y0+h-1); /* Left */ glVertex2i(x0, y0); glEnd(); glLineWidth(1.0); return; } /* End of drawRect */ /* *-------------------------------------------------------------------- * * encroached: Test a boundary edge for encroachment: the vertex * opposite the edge inside its diametral circle. * * This function is used by Function refineRd to determine which * boundary edges need to be split. * * On input: ke is the index (0 to ne-1) of a boundary edge. * * On output: The return value is 0 if ke is not a valid index of * a boundary edge. Otherwise, iv1 and iv2 index the * endpoint vertices of edge ke, and the return value * is 1 iff the vertex opposite edge ke is inside the * diametral circle; i.e., the angle opposite edge ke * in the triangle containing ke is larger than or equal & to 90 degrees. * * Global variables required by encroached: ledg, ltri, vtxy * * glmesh2 functions called by encroached: None * *-------------------------------------------------------------------- */ static int encroached(int ke, int *iv1, int *iv2) { int i, kt, kv0, kv1, kv2; if (ke < 0 || ke >= ne) return 0; kt = ledg[ke]; if (ltri[kt][6] == ke) i = 0; else if (ltri[kt][7] == ke) i = 1; else i = 2; if (ltri[kt][i+3] >= 0) return 0; kv0 = ltri[kt][i]; kv1 = ltri[kt][(i+1)%3]; kv2 = ltri[kt][(i+2)%3]; /* * Vertex kv0 encroaches upon edge ke iff <= 0, where * v0, v1, and v2 are the vertices indexed by kv0, kv1, and kv2. */ *iv1 = kv1; *iv2 = kv2; return ( (vtxy[kv1][0]-vtxy[kv0][0])* (vtxy[kv2][0]-vtxy[kv0][0]) + (vtxy[kv1][1]-vtxy[kv0][1])* (vtxy[kv2][1]-vtxy[kv0][1]) <= 0.0 ); } /* End of encroached */ /* *-------------------------------------------------------------------- * * fillRt: Fill R with triangles where R contains three adjacent * boundary vertices forming a reentrant angle (corner) at * a removable vertex. * * glmesh2 functions called: bcurve, initR, intsec * *-------------------------------------------------------------------- */ static void fillRt(void) { GLboolean intr; GLdouble x1, x2, x3, xi, y1, y2, y3, yi; GLint i, j, kc, ke, kt, kv1, kv2, kv3, nbv = 0; /* * Store listR if necessary. */ if (newR) initR(); /* * Loop on boundary curves kc. The sequence of boundary vertex * indices is stored in the first nbv elements of lste. */ for (kc = 0; kc < nc; kc++) { bcurve(kc, &nbv, lste); /* * Inner loop on three consecutive vertices (kv1,kv2,kv3) of * curve kc. */ for (j = 0; j < nbv; j++) { kv1 = lste[j]; kv2 = lste[(j+1)%nbv]; kv3 = lste[(j+2)%nbv]; /* * Test for vertices in R. */ if (!listR[kv1] || !listR[kv2] || !listR[kv3]) continue; /* * Test for kv2 removable. */ if (listn[kv2]) continue; /* * Test for kv1 adjacent to kv3. */ kt = lvtx[kv1]; if (ltri[kt][0] == kv3 || ltri[kt][1] == kv3 || ltri[kt][2] == kv3) continue; /* * Test for a reentrant corner: kv2 strictly Left kv1 -> kv3. */ x1 = vtxy[kv1][0]; y1 = vtxy[kv1][1]; x2 = vtxy[kv2][0]; y2 = vtxy[kv2][1]; x3 = vtxy[kv3][0]; y3 = vtxy[kv3][1]; if (!( (x3-x1)*(y2-y1) > (y3-y1)*(x2-x1) )) continue; /* * Test for no intersection of kv1-kv3 with an edge of the same curve * (that does not contain kv1, kv2, or kv3). */ intr = 0; i = (j+3)%nbv; while (i != (j-1+nbv)%nbv && nbv > 3) { if (intsec(x1,y1,x3,y3,vtxy[lste[i]][0],vtxy[lste[i]][1], vtxy[lste[(i+1)%nbv]][0], vtxy[lste[(i+1)%nbv]][1])) { intr = 1; break; } i = (i+1)%nbv; } if (intr) continue; if (nbv > 3) { /* * A final test for lste[i] not contained in triangle (kv1,kv3,kv2) * is necessary. */ xi = vtxy[lste[i]][0]; yi = vtxy[lste[i]][1]; if (( (x3-x1)*(yi-y1) >= (y3-y1)*(xi-x1) ) && ( (x2-x3)*(yi-y3) >= (y2-y3)*(xi-x3) ) && ( (x1-x2)*(yi-y2) >= (y1-y2)*(xi-x2) )) continue; } /* * Add triangle (kv1,kv3,kv2) and edge kv1-kv3 (unless nbv = 3). * This requires updating ltri, ledg, lvtx, nt, ne, nb, lcrv, and * lste. */ ltri[nt][0] = kv1; ltri[nt][1] = kv3; ltri[nt][2] = kv2; kt = lvtx[kv2]; if (ltri[kt][0] == kv2) i = 2; else if (ltri[kt][1] == kv2) i = 0; else i = 1; ke = ltri[kt][i+6]; ledg[ke] = nt; ltri[kt][i+3] = nt; ltri[nt][3] = kt; ltri[nt][6] = ke; kt = lvtx[kv1]; if (ltri[kt][0] == kv1) i = 2; else if (ltri[kt][1] == kv1) i = 0; else i = 1; ke = ltri[kt][i+6]; ledg[ke] = nt; ltri[kt][i+3] = nt; ltri[nt][4] = kt; ltri[nt][7] = ke; if (nbv == 3) { kt = lvtx[kv3]; if (ltri[kt][0] == kv3) i = 2; else if (ltri[kt][1] == kv3) i = 0; else i = 1; ke = ltri[kt][i+6]; ledg[ke] = nt; ltri[kt][i+3] = nt; ltri[nt][5] = kt; ltri[nt][8] = ke; nt++; nb -= 3; for (i = kc; i < nc-1; i++) { lcrv[i] = lcrv[i+1]; } nc--; break; } else { ltri[nt][5] = -1; ltri[nt][8] = ne; ledg[ne] = nt; lvtx[kv1] = nt; ne++; nt++; nb--; if (lcrv[kc] == kv2) lcrv[kc] = kv1; for (i = j+1; i < nbv-1; i++) { lste[i] = lste[i+1]; } nbv--; j--; } } } return; } /* End of fillRt */ /* *-------------------------------------------------------------------- * * filterDown: Filter down an insertion point in a binary heap in * accordance with the heap-order property: used by * buildQt and deleteQt. * * On input: ip is the index for Qtp (1 to ntq) at which the * filtering begins. * * Global variables required: ntq, Qtp * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void filterDown(int ip) { int ipn; struct triangle *p0 = Qtp[ip]; double emin0 = p0->emin; /* * Loop on tree levels. The children of Qtp[ip] are in Qtp[2*ip] * and Qtp[2*ip+1]. */ while (2*ip <= ntq) { ipn = 2*ip; if (ipn != ntq && Qtp[ipn+1]->emin < Qtp[ipn]->emin) ipn++; if (Qtp[ipn]->emin < emin0) { Qtp[ip] = Qtp[ipn]; Qtp[ip]->index_qtp = ip; } else break; ip = ipn; } Qtp[ip] = p0; Qtp[ip]->index_qtp = ip; return; } /* End of filterDown */ /* *-------------------------------------------------------------------- * * getCenter: Mouse callback (callback triggered by a mouse button * press or release with the mouse position in the * primary window) used to move the center of the viewing * volume. * * Set (xc,yc) to the object coordinates corresponding to the mouse * position when the left button is pressed. * * glmesh2 functions required: mouse, reshape, xlate * *-------------------------------------------------------------------- */ static void getCenter(int button, int state, int x, int y) { GLdouble wx, wy, wz; /* Computed object coordinates */ if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) { xlate(x,y,0.0, &wx,&wy,&wz); xc = wx; yc = wy; mode = DEFAULT; glutMouseFunc(mouse); /* Restore standard callback */ glutSetCursor(GLUT_CURSOR_INHERIT); /* Restore cursor */ reshape(glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT)); glutPostRedisplay(); } return; } /* End of getCenter */ /* *-------------------------------------------------------------------- * * getCircle: Mouse callback (callback triggered by a mouse button * press or release with the mouse position in the * primary window) used to select a circle center and * radius. * * On left button press, set cr = 0, set (cx,cy) to the object * coordinates corresponding to the mouse position, set XOR write * mode, and set mode = SELECT_CIRCLE so that Function motion allows * the user to drag a point on the circumference. * * On left button release, reset mode = DEFAULT, and restore the * standard mouse callback, inherited cursor, and overwrite mode. * * glmesh2 functions required: mouse, xlate * *-------------------------------------------------------------------- */ static void getCircle(int button, int state, int x, int y) { GLdouble wx, wy, wz; /* Computed object coordinates */ if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) { xlate(x,y,0.0, &wx,&wy,&wz); /* * Initialize circle, and set up for dragging in Function motion. */ cr = 0.0; cx = wx; cy = wy; glEnable(GL_COLOR_LOGIC_OP); glLogicOp(GL_XOR); mode = SELECT_CIRCLE; } if (button == GLUT_LEFT_BUTTON && state == GLUT_UP) { /* * Left button release: reset mode = DEFAULT, and restore the * default mouse callback, inherited cursor, and overwrite mode. */ mode = DEFAULT; glutMouseFunc(mouse); /* Restore default callback */ glutSetCursor(GLUT_CURSOR_INHERIT); /* Restore cursor */ glDisable(GL_COLOR_LOGIC_OP); glutPostRedisplay(); } return; } /* End of getCircle */ /* *-------------------------------------------------------------------- * * getCorners: Replace the set of non-removable vertices lnrv with * the set of boundary vertices that form corners (do * not lie on the line defined by their predecessors * and successors). * * glmesh2 functions called: bcurve, plDistance * *-------------------------------------------------------------------- */ static void getCorners(void) { GLint j, kc, kv, kv1, kv2, kv3, nbv = 0; GLdouble tols; /* * Compute a tolerance (squared distance) for point-line distance: * chosen to make the relative distance greater than 1.e-7. */ tols = 1.e-7*dx; tols = tols*tols; /* * Initialize nn and listn, and loop on boundary curves kc. The * sequence of boundary vertex indices is stored in the first nbv * elements of lste. */ nn = 0; for (kv = 0; kv < nv; kv++) listn[kv] = 0; for (kc = 0; kc < nc; kc++) { bcurve(kc, &nbv, lste); for (j = 0; j < nbv; j++) { kv1 = lste[(j-1+nbv)%nbv]; kv2 = lste[j]; kv3 = lste[(j+1)%nbv]; /* * For each subsequence (kv1,kv2,kv3), insert kv2 into lnrv only if * its squared distance from the line kv1-kv3 exceeds tols. */ if ( !(plDistance(vtxy[kv2][0], vtxy[kv2][1], vtxy[kv1][0], vtxy[kv1][1], vtxy[kv3][0], vtxy[kv3][1], tols) )) { lnrv[nn] = kv2; nn++; listn[kv2] = 1; } } } return; } /* End of getCorners */ /* *-------------------------------------------------------------------- * * getLin: Read a line from stdin into s. * * On input: s is a char array of length at least len. Input is * terminated by end of file, a newline character, or when * len-1 characters have been read. * * On output: s contains a (null-terminated) string of length L for * return value L <= len. Unlike the C libaray function * fgets, the newline character is not included. No * characters are left in the input stream. * * This function uses the constant EOF defined in . * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static int getLin(char *s, int len) { int c, i; i = 0; c = getchar(); while (i < len-1 && c != EOF && c != '\n') { s[i] = c; i++; c = getchar(); } s[i] = '\0'; while (c != EOF && c != '\n') { c = getchar(); /* flush input stream */ } return i+1; } /* End of getLin */ /* *-------------------------------------------------------------------- * * getPolygon: Mouse callback (callback triggered by a mouse button * press or release with the mouse position in the * primary window) used to select a polygonal region R. * * On first left button press (np = 0), set (xv1,yv1) and (xv2,yv2) * to the object coordinates corresponding to the mouse position when * the left button is pressed, store (xv1,yv1) in R, and set XOR * write mode. On input, mode = CONSTRUCT_R so that the mouse motion * callback (Function motion) drags a line segment with the mouse * position. * * On subsequent left button presses, set (xv2,yv2) to the mouse * position, and draw the line segment from (xv1,yv1) to (xv2,yv2). * * On left button release, set (xv2,yv2) to the mouse position, add * the line segment from (xv1,yv1) to (xv2,yv2) to R, and copy * (xv2,yv2) to (xv1,yv1) unless the new line segment would result in * a self-intersecting curve. * * glmesh2 functions required: intsec, xlate * *-------------------------------------------------------------------- */ static void getPolygon(int button, int state, int x, int y) { GLdouble wx, wy, wz; /* Computed world coordinates */ GLint i; GLboolean intr; if (button != GLUT_LEFT_BUTTON) return; /* * Left button: compute object coordinates (wx,wy) corresponding * to the mouse position. */ xlate(x,y,0.0, &wx,&wy,&wz); if (state == GLUT_DOWN && np == 0) { /* * First left button press: store (wx,wy) in (xv1,yv1), (xv2,yv2), * and R. */ xv1 = wx; yv1 = wy; xv2 = wx; yv2 = wy; R[np][0] = wx; R[np][1] = wy; np++; /* * Set XOR write mode. */ glEnable(GL_COLOR_LOGIC_OP); glLogicOp(GL_XOR); } if (state == GLUT_DOWN && np > 0) { /* * Left button press: set (xv2,yv2) to (wx,wy), and * draw line segment from (xv1,xv2) to (xv2,yv2). */ xv2 = wx; yv2 = wy; glBegin(GL_LINES); glVertex2d(xv1,yv1); glVertex2d(xv2,yv2); glEnd(); glutSwapBuffers(); } if (state == GLUT_UP) { /* * Left button release: set (xv2,yv2) to (wx,wy). */ xv2 = wx; yv2 = wy; /* * Test for an intersection. */ if (xv2 == xv1 && yv2 == yv1) return; intr = 0; for (i = 0; i <= np-3; i++) { if (intsec(xv1,yv1,xv2,yv2,R[i][0],R[i][1],R[i+1][0], R[i+1][1])) { intr = 1; break; } } if (intr) { /* * Invalid choice: erase the line segment, and reset (xv2,yv2) to * (xv1,yv1). */ glBegin(GL_LINES); glVertex2d(xv1,yv1); glVertex2d(xv2,yv2); glEnd(); glutSwapBuffers(); xv2 = xv1; yv2 = yv1; return; } /* * Store (xv2,yv2) in R, and copy (xv2,yv2) to (xv1,yv1). */ if (np >= npmax) { npmax *= 2; R = realloc(R, npmax*sizeof *R); if (R == NULL) { printf("glmesh2: Unable to allocate sufficient memory" " for R.\n"); exit(1); } } R[np][0] = xv2; R[np][1] = yv2; np++; xv1 = xv2; yv1 = yv2; } return; } /* End of getPolygon */ /* *-------------------------------------------------------------------- * * getTitlePos: Mouse callback (callback triggered by a mouse button * press or release with the mouse position in the * primary window) used to select a title position. * * This function sets title_pos to the object coordinates of the * mouse position when the left button is pressed. * * glmesh2 functions required: mouse, xlate * *-------------------------------------------------------------------- */ static void getTitlePos(int button, int state, int x, int y) { GLdouble wx, wy, wz; /* Computed object coordinates */ if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) { xlate(x,y,0.0, &wx,&wy,&wz); title_pos[0] = (GLfloat) wx; title_pos[1] = (GLfloat) wy; title_pos[2] = (GLfloat) wz; mode = DEFAULT; glutMouseFunc(mouse); /* Restore standard callback */ glutSetCursor(GLUT_CURSOR_INHERIT); /* Restore cursor */ glutPostRedisplay(); } return; } /* End of getTitlePos */ /* *-------------------------------------------------------------------- * * getVertex: Mouse callback (callback triggered by a mouse button * press or release with the mouse position in the * primary window) used to select a vertex in uv-edit * mode. * * On left button press, if the mouse position is close enough to a * vertex, write the vertex index, x, y, u, and v values to stdout, * and prompt for new u and v values. * * glmesh2 function called: xlate * *-------------------------------------------------------------------- */ static void getVertex(int button, int state, int x, int y) { GLdouble wx, wy, wz; /* Computed world coordinates */ GLboolean error, intr; GLdouble tols; GLint kv; int linsiz = 80; char line[80]; char *readerror = "Invalid entry.\n"; if (button != GLUT_LEFT_BUTTON || state != GLUT_DOWN) return; /* * Left button press: compute the object coordinates (wx,wy) of the * mouse position. */ xlate(x,y,0.0, &wx,&wy,&wz); /* * Find a vertex kv, if any, within squared distance tols of the * mouse position. */ tols = (sx <= sy) ? sx : sy; /* tols = min(sx, sy) */ tols = 25.0*tols*tols; /* 5-pixel tolerance */ intr = 0; for (kv = 0; kv < nv; kv++) { if ((intr = ((wx-vtxy[kv][0])*(wx-vtxy[kv][0]) + (wy-vtxy[kv][1])*(wy-vtxy[kv][1]) <= tols))) break; } if (!intr) return; /* * Write kv, x, y, u, and v to stdout, and prompt for new (u,v). */ printf("\n Vertex %d: x = %e, y = %e, u = %e, v = %e\n", kv, vtxy[kv][0], vtxy[kv][1], uv[kv][0], uv[kv][1]); error = 1; while (error) { printf(" Enter u: "); if (fgets(line, linsiz, stdin) == NULL) return; error = (sscanf(line, "%le", &uv[kv][0]) != 1); if (error) printf("%s", readerror); } error = 1; while (error) { printf(" Enter v: "); if (fgets(line, linsiz, stdin) == NULL) return; error = (sscanf(line, "%le", &uv[kv][1]) != 1); if (error) printf("%s", readerror); } printf(" Vertex %d: x = %e, y = %e, u = %e, v = %e\n\n", kv, vtxy[kv][0], vtxy[kv][1], uv[kv][0], uv[kv][1]); return; } /* End of getVertex */ /* *-------------------------------------------------------------------- * * getViewVolume: Mouse callback (callback triggered by a mouse * button press or release with the mouse position * in the primary window) used to select a new * viewing volume (pair of corner coordinates). * * On left button press, set (xv1,yv1) and (xv2,yv2) to the object * coordinates corresponding to the mouse position when the left * button is pressed, set XOR write mode, and set mode = ALTER_VIEW_ * VOLUME so that the mouse motion callback (Function motion) drags * a rectangle outline with the mouse position. * * On left button release, set (xv2,yv2) to the mouse position, set * (xc,yc), dx, and dy to the center, width, and height of the * rectangle [xv1,xv2] X [yv1,yv2], restore the standard callback, * cursor, and write mode, and call reshape. * * glmesh2 functions required: mouse, reshape, xlate * *-------------------------------------------------------------------- */ static void getViewVolume(int button, int state, int x, int y) { GLdouble wx, wy, wz; /* Computed world coordinates */ if (button != GLUT_LEFT_BUTTON) return; /* * Left button: compute object coordinates (wx,wy) corresponding * to the mouse position. */ xlate(x,y,0.0, &wx,&wy,&wz); if (state == GLUT_DOWN) { /* * Left button press: set (xv1,yv1) and (xv2,yv2) to (wx,wy). */ xv1 = wx; yv1 = wy; xv2 = xv1; yv2 = yv1; /* * Set XOR write mode. */ glEnable(GL_COLOR_LOGIC_OP); glLogicOp(GL_XOR); /* * Set mode = ALTER_VIEW_VOLUME and disable the cursor. */ mode = ALTER_VIEW_VOLUME; glutSetCursor(GLUT_CURSOR_NONE); } else { /* * Left button release: reset mode to DEFAULT and set (xv2,yv2) to * (wx,wy). */ mode = DEFAULT; xv2 = wx; yv2 = wy; /* * Restore standard mouse callback, inherited cursor, and overwrite * mode. */ glutMouseFunc(mouse); glutSetCursor(GLUT_CURSOR_INHERIT); glDisable(GL_COLOR_LOGIC_OP); if (xv2 == xv1 || yv2 == yv1) return; /* * Compute xc, yc, dx, and dy, call reshape, and post redisplay. */ xc = (xv1+xv2)/2.0; yc = (yv1+yv2)/2.0; dx = xv2-xv1; if (dx < 0.0) dx = -dx; dy = yv2-yv1; if (dy < 0.0) dy = -dy; reshape(glutGet(GLUT_WINDOW_WIDTH),glutGet(GLUT_WINDOW_HEIGHT)); glutPostRedisplay(); } return; } /* End of getViewVolume */ /* *-------------------------------------------------------------------- * * init0: Initialization for window 0. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void init0(void) { /* * Set background color. */ glClearColor (color_b[0], color_b[1], color_b[2], 0.0); glShadeModel(GL_FLAT); /* Default is smooth */ return; } /* End of init0 */ /* *-------------------------------------------------------------------- * * initR: Initialization for user-selected polygon R: store flags * listR. * * Global variables required: * * nv = Numbers of vertices in the mesh. * * vtxy = Vertex coordinates. * * newR = Flag with value TRUE iff Function initR needs to be * called before the triangle mesh can be altered: set * to FALSE here. * * listR = Array of length nvmax used to store Boolean flags * marking vertices as contained in R. * * np = Number of vertices in the boundary of R. * * R = Array dimensioned npmax by 2 containing the Cartesian * coordinates of the sequence of vertices defining the * boundary of R. * * glmesh2 function called: inside * *-------------------------------------------------------------------- */ static void initR(void) { int kv; /* * Test for at least three boundary vertices. */ newR = 0; if (np < 3) { for (kv = 0; kv < nv; kv++) listR[kv] = 0; return; } /* * Mark the vertices that are found in R by storing flags in listR. */ for (kv = 0; kv < nv; kv++) { if (inside(vtxy[kv][0],vtxy[kv][1])) { listR[kv] = 1; } else { listR[kv] = 0; } } return; } /* End of initR */ /* *-------------------------------------------------------------------- * * inputData: Allocate storage and read a data set containing a tri- * angulation: number of vertices nv, vertex coordinates * vtxy, number of triangles nt, and CCW-ordered triples * of triangle vertex indices (first three columns of * ltri). * * On input: argc and argv are the parameters input to main. * * On output: Unless the program is terminated with exit code 1 * (and an error message), the global variables defining * the triangle mesh are initialized with data from the * input file. Storage is also allocated for lste, listn, * and listR, and output file name fname is stored. * * Error conditions are tested in the following order: * * 1: read error. This condition is tested following every * read. * 2: nv < 3. * 3: unable to allocate sufficient memory. This condition * is tested following every call to malloc. * 4: nt < nv-2. * 5: an ltri vertex index is outside its valid range * (0 to nv-1). * 6: a triangle's vertex indices are not distinct. * 7: two triangles have the same directed edge (pair of * adjacent vertex indices). * 8: a set of three vertex indices occurs in both clockwise * and CCW order in ltri. * 9: the same sequence of vertex indices occurs more than once * in ltri (duplicate triangles). * 10: a triangle has two boundary edges. This is a problem for * Dirichlet boundary conditions but not on an outflow * boundary. Therefore, a warning message is written, and * execution is continued. * 11: an open boundary curve was encountered. * 12: the computed parameters ne, nc, and nb are not consistent * with nv and nt. This could be caused by a vertex index * appearing more than once in the set of boundary curves. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void inputData(int argc, char *argv[]) { int i, i0, i1, j, k, ke, kt, kt0, ktn, kts; int linsiz = 160; char line[160]; FILE *fpr; char *filemsg1 = "\n\n" " glmesh2: 2-D triangle mesh generation using OpenGL.\n\n" " This program requires two command-line arguments: an input\n" " file name (or path) and an output file name, where the\n" " input file contains the following records:\n\n" " Optional comment lines beginning with # in position 1\n" " nv = number of vertices (at least 3)\n" " x y = Cartesian coordinates of vertex 0\n" " x y = Cartesian coordinates of vertex 1\n" " . . .\n" " x y = Cartesian coordinates of vertex nv-1\n"; char *filemsg2 = " Optional comment lines beginning with # in position 1\n" " nt = number of triangles (at least nv-2)\n" " i1 i2 i3 = CCW-ordered vertex indices of triangle 0\n" " i1 i2 i3 = CCW-ordered vertex indices of triangle 1\n" " . . .\n" " i1 i2 i3 = CCW-ordered vertex indices of triangle nt-1\n\n" " Vertex indices are in the range 0 to nv-1.\n"; char *readerror = "glmesh2: Error reading file.\n"; char *memerror = "glmesh2: Unable to allocate sufficient memory.\n"; /* * Open input file with pointer fpr, and store argv[2] in fname. */ if (argc < 3) { printf("%s", filemsg1); printf("%s", filemsg2); exit(1); } fpr = fopen(argv[1], "r"); if (fpr == NULL) { printf("glmesh2: Cannot open file %s.\n", argv[1]); exit(1); } fname = argv[2]; /* * Read a line at a time from fpr into line, beginning with comments. */ line[0] = '#'; while (line[0] == '#') if (fgets(line, linsiz, fpr) == NULL) { printf("%s", readerror); exit(1); } /* * Read nv. */ if (sscanf(line, "%d", &nv) != 1) { printf("%s", readerror); exit(1); } if (nv < 3) { printf("glmesh2: Invalid value: nv = %d.\n", nv); exit(1); } /* * Initialize nvmax to 4*nv, and allocate storage for vtxy, uv, lvtx, * ltri, ledg, lnrv, lcrv, lste, listn, and listR. */ nvmax = 4*nv; vtxy = malloc(nvmax*sizeof *vtxy); uv = malloc(nvmax*sizeof *uv); lvtx = malloc(nvmax*sizeof *lvtx); ltri = malloc(2*nvmax*sizeof *ltri); ledg = malloc(3*nvmax*sizeof *ledg); lnrv = malloc(nvmax*sizeof *lnrv); lcrv = malloc((nvmax/3)*sizeof *lcrv); lste = malloc(3*nvmax*sizeof *lste); listn = malloc(nvmax*sizeof *listn); listR = malloc(nvmax*sizeof *listR); if (vtxy == NULL || uv == NULL || lvtx == NULL || ltri == NULL || ledg == NULL || lnrv == NULL || lcrv == NULL || lste == NULL || listn == NULL || listR == NULL) { printf("%s", memerror); exit(1); } /* * Read vertices vtxy. */ for (i = 0; i < nv; i++) { if (fgets(line, linsiz, fpr) == NULL) { printf("%s", readerror); exit(1); } if (sscanf(line, "%lf%lf", &vtxy[i][0], &vtxy[i][1]) != 2) { printf("%s", readerror); exit(1); } } /* * Store zeros in uv. */ for (i = 0; i < nv; i++) { uv[i][0] = 0.0; uv[i][1] = 0.0; } /* * Read comment lines if any. */ line[0] = '#'; while (line[0] == '#') if (fgets(line, linsiz, fpr) == NULL) { printf("%s", readerror); exit(1); } /* * Read number of triangles nt. */ if (sscanf(line, "%d", &nt) != 1) { printf("%s", readerror); exit(1); } if (nt < nv-2) { printf("glmesh2: Invalid value: nt = %d.\n", nt); exit(1); } /* * Read triangle list ltri. */ for (kt = 0; kt < nt; kt++) { if (fgets(line, linsiz, fpr) == NULL) { printf("%s", readerror); exit(1); } if (sscanf(line, "%d%d%d", <ri[kt][0], <ri[kt][1], <ri[kt][2]) != 3) { printf("%s", readerror); exit(1); } if (ltri[kt][0] < 0 || ltri[kt][0] >= nv || ltri[kt][1] < 0 || ltri[kt][1] >= nv || ltri[kt][2] < 0 || ltri[kt][2] >= nv) { printf("glmesh2: The i-th triple includes an invalid index " "for i = %d.\n", kt); exit(1); } if (ltri[kt][0] == ltri[kt][1] || ltri[kt][1] == ltri[kt][2] || ltri[kt][2] == ltri[kt][0]) { printf("glmesh2: The vertex indices are not distinct in the " "i-th triple for i = %d.\n", kt); exit(1); } } /* * Initialize ltri columns 3 to 8 with flag values -1 to indicate * values that have not yet been computed, and initialize counts. */ for (kt = 0; kt < nt; kt++) { for (j = 3; j < 9; j++) { ltri[kt][j] = -1; } } nc = 0; nb = 0; ne = 0; /* * Construct the remaining portion of ltri by searching for * neighboring triangles and generating edge indices. The outer * loop is on triangles kt; inner loop on potential neighbors ktn. */ for (kt = 1; kt < nt; kt++) { for (ktn = 0; ktn < kt; ktn++) { /* * Loop on pairs of edges ltri[kt][i]-ltri[kt][i+1], * ltri[ktn][j]-ltri[ktn][j+1]. */ for (i = 0; i < 3; i++) { for (j = 0; j < 3; j++) { /* * Test for error 7. */ if (ltri[kt][i] == ltri[ktn][j] && ltri[kt][(i+1)%3] == ltri[ktn][(j+1)%3]) { printf("glmesh2: Two or more triangles have the " "same pair of adjacent vertex indices.\n"); exit(1); } if (ltri[kt][i] != ltri[ktn][(j+1)%3] || ltri[kt][(i+1)%3] != ltri[ktn][j]) continue; /* * Shared edge: ltri[kt][i]-ltri[kt][i+1] = * ltri[ktn][j+1]-ltri[ktn][j] * * Test for errors 7, 8, and 9, and store the indices. */ if (ltri[kt][(i+2)%3] == ltri[ktn][(j+2)%3]) { printf("glmesh2: The same triple of vertex indices " "occurs\n with both cyclical orderings.\n"); exit(1); } if (ltri[ktn][(j+2)%3+3] >= 0) { k = ltri[ktn][(j+2)%3+3]; if (ltri[kt][(i+2)%3] == ltri[k][(i+2)%3]) { printf("glmesh2: Duplicate vertex index " "triples.\n"); exit(1); } printf("glmesh2: Two or more triangles have the " "same pair of adjacent vertex indices.\n"); exit(1); } ltri[kt][(i+2)%3+3] = ktn; ltri[ktn][(j+2)%3+3] = kt; ne++; ltri[kt][(i+2)%3+6] = ne-1; ltri[ktn][(j+2)%3+6] = ne-1; } } } } /* * Count the boundary curves nc and boundary edges (vertices) nb, * assign them edge indices, and construct lcrv. */ for (kt0 = 0; kt0 < nt; kt0++) { for (j = 0; j < 3; j++) { if (ltri[kt0][j+6] >= 0) continue; /* * New boundary curve beginning with vertices i0, i1. * Initialize for traversal. A counter k is used to * avoid an infinite loop. */ kts = kt0; kt = kt0; i = j; nc++; nb++; ne++; ltri[kt][i+6] = ne-1; i = (i+1)%3; i1 = ltri[kt][i]; i = (i+1)%3; i0 = ltri[kt][i]; lcrv[nc-1] = i0; k = 1; /* * Boundary curve traversal: find the next vertex. */ while (i1 != i0) { ktn = ltri[kt][i+3]; while (ktn >= 0) { kt = ktn; if (ltri[kt][0] == i1) i = 1; else if (ltri[kt][1] == i1) i = 2; else i = 0; ktn = ltri[kt][i+3]; } if (kt == kts) printf("glmesh2: WARNING: Triangle %d has two " "boundary edges.\n", kt); kts = kt; nb++; ne++; ltri[kt][i+6] = ne-1; i = (i+1)%3; i1 = ltri[kt][i]; i = (i+1)%3; k++; if (k > nv) { printf("glmesh2: Open boundary curve.\n"); exit(1); } } } } /* * Test for error 12. */ if (nb < 3*nc || ne != 3*nv-nb+3*nc-6 || nt != 2*nv-nb+2*nc-4) { printf("glmesh2: Inconsistent counts: nv = %d, nb = %d, " "nc = %d, nt = %d, ne = %d\n", nv, nb, nc, nt, ne); exit(1); } /* * Construct lvtx. Initialize the list with flag values -1 * indicating entries not yet set. */ for (i0 = 0; i0 < nv; i0++) lvtx[i0] = -1; /* * Store triangle indices of boundary vertices i0 in a loop * on triangle/vertex pairs. */ for (kt = 0; kt < nt; kt++) { for (i = 0; i < 3; i++) { if (ltri[kt][i+3] < 0) { i0 = ltri[kt][(i+1)%3]; lvtx[i0] = kt; } } } /* * Fill in triangle indices for vertices with negative entries, * and store ledg. */ for (kt = 0; kt < nt; kt++) { for (i = 0; i < 3; i++) { i0 = ltri[kt][i]; if (lvtx[i0] < 0) lvtx[i0] = kt; if (ltri[kt][i+3] < kt) { ke = ltri[kt][i+6]; ledg[ke] = kt; } } } return; } /* End of inputData */ /* *-------------------------------------------------------------------- * * insertQt: Insert a triangle into priority queue Qt. * * On input: kt is the index (0 to nt-1) of the triangle to be * inserted, imin is the local index (0, 1, or 2) of the * vertex opposite the shortest edge of triangle kt, and * emin is the squared length of the shortest edge (the * key value). * * Global variables required: ltri, ntq, ntqmax, Qt, Qtp * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void insertQt(int kt, int imin, double emin) { int i, ip; struct triangle *pnew; /* * Test for sufficient storage in Qt. Note that, since the address * of Qt is likely to be changed by realloc, the pointers in Qtp * must be restored. */ if (ntq == ntqmax) { ntqmax *= 2; Qt = realloc(Qt, ntqmax*sizeof *Qt); Qtp = realloc(Qtp, (ntqmax+1)*sizeof *Qtp); for (i = 0; i < ntq; i++) Qtp[Qt[i].index_qtp] = &Qt[i]; } /* * Store a new triangle in Qt[ntq] with pointer pnew, and * increment ntq. */ Qt[ntq].index_qtp = ntq+1; Qt[ntq].index_t = kt; Qt[ntq].index_v1 = ltri[kt][0]; Qt[ntq].index_v2 = ltri[kt][1]; Qt[ntq].index_v3 = ltri[kt][2]; Qt[ntq].imin = imin; Qt[ntq].emin = emin; pnew = &Qt[ntq]; ntq++; /* * The insertion point ip, initially the first empty array position, * is percolated up as required by the heap-order property. */ ip = ntq; while (ip > 1 && emin < Qtp[ip/2]->emin) { Qtp[ip] = Qtp[ip/2]; Qtp[ip]->index_qtp = ip; ip /= 2; } Qtp[ip] = pnew; Qtp[ip]->index_qtp = ip; return; } /* End of insertQt */ /* *-------------------------------------------------------------------- * * insertv: Partition a triangle into three subtriangles by insert- * ing a new vertex in its interior. The uv values at the * new vertex are computed by linear interpolation. * * On input: kt0 is the index of the triangle to be partitioned. * * b1, b2, and b3 are the barycentric coordinates of the * new vertex with respect to triangle kt0. It is assumed * without a test that these values are nonnegative and * add to 1. * * On output: Unless nv >= nvmax or kt0 is an invalid index, the * triangle mesh is updated with the revision. The new * vertex is flagged as being contained in R (listR[nv-1] * = 1) iff all three vertices of triangle kt0 are in R. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void insertv(int kt0, double b1, double b2, double b3) { int i, ke1, ke2, ke3, kt1, kt2, kt3, kv1, kv2, kv3; /* * Test for invalid input. */ if (nv >= nvmax || kt0 < 0 || kt0 >= nt) return; /* * Store indices of vertices, opposite triangles, and edges of kt0 * in local variables. */ kv1 = ltri[kt0][0]; kv2 = ltri[kt0][1]; kv3 = ltri[kt0][2]; kt1 = ltri[kt0][3]; kt2 = ltri[kt0][4]; kt3 = ltri[kt0][5]; ke1 = ltri[kt0][6]; ke2 = ltri[kt0][7]; ke3 = ltri[kt0][8]; /* * Compute the Cartesian coordinates of the new vertex (index nv), * and update vtxy and uv. */ vtxy[nv][0] = b1*vtxy[kv1][0] + b2*vtxy[kv2][0] + b3*vtxy[kv3][0]; vtxy[nv][1] = b1*vtxy[kv1][1] + b2*vtxy[kv2][1] + b3*vtxy[kv3][1]; uv[nv][0] = b1*uv[kv1][0] + b2*uv[kv2][0] + b3*uv[kv3][0]; uv[nv][1] = b1*uv[kv1][1] + b2*uv[kv2][1] + b3*uv[kv3][1]; /* * Overwrite row kt0, append two new rows to ltri, and adjust * additional ltri entries. */ ltri[kt0][0] = nv; ltri[kt0][4] = nt; ltri[kt0][5] = nt+1; ltri[kt0][7] = ne+2; ltri[kt0][8] = ne+1; ltri[nt][0] = nv; ltri[nt][1] = kv3; ltri[nt][2] = kv1; ltri[nt][3] = kt2; ltri[nt][4] = nt+1; ltri[nt][5] = kt0; ltri[nt][6] = ke2; ltri[nt][7] = ne; ltri[nt][8] = ne+2; ltri[nt+1][0] = nv; ltri[nt+1][1] = kv1; ltri[nt+1][2] = kv2; ltri[nt+1][3] = kt3; ltri[nt+1][4] = kt0; ltri[nt+1][5] = nt; ltri[nt+1][6] = ke3; ltri[nt+1][7] = ne+1; ltri[nt+1][8] = ne; if (kt2 >= 0) { if (ltri[kt2][3] == kt0) i = 3; else if (ltri[kt2][4] == kt0) i = 4; else i = 5; ltri[kt2][i] = nt; } if (kt3 >= 0) { if (ltri[kt3][3] == kt0) i = 3; else if (ltri[kt3][4] == kt0) i = 4; else i = 5; ltri[kt3][i] = nt+1; } /* * Adjust ledg and lvtx entries. */ ledg[ke2] = nt; ledg[ke3] = nt+1; ledg[ne] = nt+1; ledg[ne+1] = nt+1; ledg[ne+2] = nt; if (lvtx[kv1] == kt0) lvtx[kv1] = nt+1; if (lvtx[kv3] == kt0) lvtx[kv3] = nt; lvtx[nv] = kt0; /* * Store listn[nv] and listR[nv]. */ listn[nv] = 0; listR[nv] = listR[kv1] && listR[kv2] && listR[kv3]; /* * Update vertex, edge, and triangle counts. */ nv++; ne += 3; nt += 2; return; } /* End of insertv */ /* *-------------------------------------------------------------------- * * inside: Locate a point P relative to the polygon R, returning a * nonzero value iff P is contained in R. * * R is defined by a cyclically ordered sequence of vertices defining * a simple closed curve. Adjacent vertices need not be distinct, * and the curve orientation is irrelevant. * * The algorithm is referred to as ray casting, crossing number, or * even-odd rule. The point P is inside R iff there are an odd * number of intersections between the polygon boundary curve and a * horizontal ray with endpoint P. Intersections are counted by * testing each curve segment (polygon side) for intersection. In * order to avoid counting a vertex as two intersections when it lies * on the ray and has neighbors on opposite sides of the ray, we * arbitrarily define vertices with the same y component as P as * being above the ray. * * If P lies on the boundary of R, the decision is unspecified, and * if P is close to the polygon boundary, the problem is ill- * conditioned, and the decision may be incorrect. * * On input: xp,yp = Cartesian coordinates of P. * * Global variables required: * * np = Number of vertices in the boundary of R. * * R = Array dimensioned np by 2 containing the Cartesian * coordinates of the sequence of vertices defining the * boundary of R. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static unsigned char inside(double xp, double yp) { unsigned char odd = 0; int i; double x, x1, x2, y1, y2; x2 = R[np-1][0]; y2 = R[np-1][1]; for (i = 0; i < np; i++) { x1 = x2; y1 = y2; x2 = R[i][0]; y2 = R[i][1]; /* * Test for an intersection with the line y = yp. */ if ((y1 < yp && y2 < yp) || (y1 >= yp && y2 >= yp)) continue; /* * Compute the x component of the point of intersection, and test * for an intersection with the ray x > xp. */ x = x1 + (yp-y1)*(x2-x1)/(y2-y1); if (x > xp) odd = !odd; } return odd; } /* End of inside */ /* *-------------------------------------------------------------------- * * intsec: Test for an intersection between a pair of line * segments in the plane. * * Note that an incorrect decision may result from floating point * error if the four endpoints are nearly collinear. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static unsigned char intsec(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4) { double a, b, d, dx12, dx31, dx34, dy12, dy31, dy34; /* * Test for overlap between the smallest rectangles that contain * the line segments and have sides parallel to the axes. */ if ((x1 < x3 && x1 < x4 && x2 < x3 && x2 < x4) || (x1 > x3 && x1 > x4 && x2 > x3 && x2 > x4) || (y1 < y3 && y1 < y4 && y2 < y3 && y2 < y4) || (y1 > y3 && y1 > y4 && y2 > y3 && y2 > y4)) { return 0; } /* * Compute a = (p4-p3) X (p1-p3), b = (p2-p1) X (p1-p3), and * d = (p2-p1) X (p4-p3) (z components). */ dx12 = x2 - x1; dy12 = y2 - y1; dx34 = x4 - x3; dy34 = y4 - y3; dx31 = x1 - x3; dy31 = y1 - y3; a = dx34*dy31 - dx31*dy34; b = dx12*dy31 - dx31*dy12; d = dx12*dy34 - dx34*dy12; /* * If d == 0, either the line segments are parallel, or one (or both) * of them is a single point. */ if (d == 0.0) return (a == 0.0 && b == 0.0); /* * d != 0 and the point of intersection of the lines defined by the * line segments is p = p1 + (a/d)*(p2-p1) = p3 + (b/d)*(p4-p3). */ return (a/d >= 0.0 && a/d <= 1.0 && b/d >= 0.0 && b/d <= 1.0); } /* End of intsec */ /* *-------------------------------------------------------------------- * * key: Associate ASCII keys with menu entries. * The mouse coordinates (x,y) are not used. * * glmesh2 function called: menu * *-------------------------------------------------------------------- */ static void key(unsigned char key, int x, int y) { menu((int) key); return; } /* End of key */ /* *-------------------------------------------------------------------- * * makeMenu: Create a menu (for window 0) attached to the right * mouse button. * * glmesh2 function required: menu * *-------------------------------------------------------------------- */ static void makeMenu(void) { int submenu1, submenu2, submenu3; glutCreateMenu(menu); glutAddMenuEntry("", 0); glutAddMenuEntry(" B: Set R to contain the Bounding box", 'B'); glutAddMenuEntry(" c: Cycle through triangle Color-fill options", 'c'); glutAddMenuEntry(" e: Toggle annotation of Edges with indices", 'e'); glutAddMenuEntry(" F: Output the triangle mesh to a File", 'F'); glutAddMenuEntry(" o: Draw a circle (drag with left mouse button)", 'o'); glutAddMenuEntry(" P: Print triangle mesh data structure", 'P'); glutAddMenuEntry(" q: Query current state (print parameters)", 'q'); glutAddMenuEntry(" r: Restore default viewing volume", 'r'); glutAddMenuEntry(" R: Select a new polygon R (mouse button presses)", 'R'); glutAddMenuEntry(" t: Toggle annotation of Triangles with indices", 't'); glutAddMenuEntry(" u: Enter uv-edit mode (select vertex with mouse)", 'u'); glutAddMenuEntry(" v: Toggle annotation of Vertices with indices", 'v'); glutAddMenuEntry(" V: Toggle annotation of non-removable Vertices", 'V'); glutAddMenuEntry(" z: Zoom in", 'z'); glutAddMenuEntry(" Z: Zoom out", 'Z'); glutAddMenuEntry(" [Arrows]: Move center", 0); glutAddMenuEntry(" [Left mouse button]: Alter polygon R", 0); glutAddMenuEntry(" [Shifted left mouse button]: Delete vertex of R", 0); glutAddMenuEntry(" [Enter]: Terminate mode (D, E, or R)", 13); glutAddMenuEntry(" [esc]: Cancel mode (D/E/R), or terminate program", 27); submenu1 = glutCreateMenu(menu); glutAddMenuEntry(" a: Toggle preservation of Aspect ratio", 'a'); glutAddMenuEntry(" b: Toggle display of Bounding box", 'b'); glutAddMenuEntry(" H: Toggle off-center/circumcenter in refineRd (&)",'H'); glutAddMenuEntry(" i: Triangulation Integrity check (debugging)", 'i'); glutAddMenuEntry(" j: Toggle antialiasing (Jaggies)", 'j'); glutAddMenuEntry(" k: Decrease line width", 'k'); glutAddMenuEntry(" K: Increase line width", 'K'); glutAddMenuEntry(" m: Decrease minimum angle for refineRd (&)", 'm'); glutAddMenuEntry(" M: Increase minimum angle for refineRd (&)", 'M'); glutAddMenuEntry(" p: Pan (select new center with mouse button)", 'p'); glutAddMenuEntry(" w: Select a new view volume (left mouse button)", 'w'); glutAddMenuEntry(" x: Halve max. no. Steiner pts. for refineRd (&)", 'x'); glutAddMenuEntry(" X: Double max. no. Steiner pts. for refineRd (&)", 'X'); glutAddMenuEntry(" [Home]: Restore default viewing volume", 0); glutAddMenuEntry(" [Page Up, Page Down]: Zoom in, zoom out", 0); submenu2 = glutCreateMenu(menu); glutAddMenuEntry(" C: Replace non-removable vertex set by Corners", 'C'); glutAddMenuEntry(" D: Alter Domain (drag non-removable vertices)", 'D'); glutAddMenuEntry(" E: Select Edges for swapping (left mouse button)", 'E'); glutAddMenuEntry(" f: Fill concavities with triangles", 'f'); glutAddMenuEntry(" L: Toggle boundary Layer (smoothing method)", 'L'); glutAddMenuEntry(" n: Convert vertices in R to removable", 'n'); glutAddMenuEntry(" N: Convert vertices in R to Non-removable", 'N'); glutAddMenuEntry(" O: Move non-removable vertices in R to circle", 'O'); glutAddMenuEntry(" s: Optimize R with Delaunay edge Swaps", 's'); glutAddMenuEntry(" S: Optimize R by Smoothing (moving vertices)", 'S'); glutAddMenuEntry(" &: Refine R using Delaunay refinement", '&'); glutAddMenuEntry(" +: Refine R by adding triangle barycenters", '+'); glutAddMenuEntry(" >: Refine R by adding edge midpoints", '>'); glutAddMenuEntry(" -: Unrefine by reversing &, +, or > if possible", '-'); glutAddMenuEntry(" <: Coarsen R by removing vertices", '<'); glutAddMenuEntry(" [Backspace or Delete]: Remove triangles in R", 8); submenu3 = glutCreateMenu(menu); glutAddMenuEntry(" T: Replace Text (keyboard string entry)", 'T'); glutAddMenuEntry(" l: Change title Location (left button)", 'l'); glutAddMenuEntry(" !: Toggle font selection: indices/title", '!'); glutAddMenuEntry(" 0: Set 8 by 13 fixed width font", '0'); glutAddMenuEntry(" 1: Set 9 by 15 fixed width font", '1'); glutAddMenuEntry(" 2: Set 10-point prop. spaced Times Roman", '2'); glutAddMenuEntry(" 3: Set 24-point prop. spaced Times Roman", '3'); glutAddMenuEntry(" 4: Set 10-point prop. spaced Helvetica", '4'); glutAddMenuEntry(" 5: Set 12-point prop. spaced Helvetica", '5'); glutAddMenuEntry(" 6: Set 18-point prop. spaced Helvetica", '6'); glutAddMenuEntry(" 7: Set prop. spaced Roman stroke font", '7'); glutAddMenuEntry(" 8: Set mono-spaced Roman stroke font", '8'); glutSetMenu(1); glutAddSubMenu(" Additional Options -->", submenu1); glutAddSubMenu(" Triangulation Updates -->", submenu2); glutAddSubMenu(" Title and Font -->", submenu3); glutAttachMenu(GLUT_RIGHT_BUTTON); return; } /* End of makeMenu */ /* *-------------------------------------------------------------------- * * menu: Respond to menu selections. * * glmesh2 functions required: all of them * *-------------------------------------------------------------------- */ static void menu(int item) { double c; static unsigned char font_select = 0; switch (item) { case 'a': /* Toggle preservation of aspect ratio */ aspect = !aspect; reshape(glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT)); break; case 'b': /* Toggle display of bounding box */ plot_box = !plot_box; break; case 'B': /* Set R to a rectangle containing the bounding box */ np = 4; R[0][0] = xmin - 0.05*(xmax-xmin); R[0][1] = ymin - 0.05*(ymax-ymin); R[1][0] = xmax + 0.05*(xmax-xmin); R[1][1] = ymin - 0.05*(ymax-ymin); R[2][0] = xmax + 0.05*(xmax-xmin); R[2][1] = ymax + 0.05*(ymax-ymin); R[3][0] = xmin - 0.05*(xmax-xmin); R[3][1] = ymax + 0.05*(ymax-ymin); initR(); break; case 'c': /* Increment color-fill flag */ color_fill = (color_fill+1)%4; break; case 'C': /* Replace set of non-removable vertices lnrv by the set of corners */ getCorners(); break; case 'D': /* Drag a non-removable vertex */ if (mode == CONSTRUCT_R) { /* Terminate polygon selection */ np = 0; } if (newR) initR(); mode = ALTER_DOMAIN; glutMouseFunc(dragCorner); glutSetCursor(GLUT_CURSOR_CROSSHAIR); /* Change cursor */ break; case 'e': /* Toggle display of edge indices */ plot_e = !plot_e; break; case 'E': /* Enter edge-swapping mode (left button selection) */ if (mode == CONSTRUCT_R) { /* Terminate polygon selection */ glDisable(GL_COLOR_LOGIC_OP); np = 0; } if (mode == ALTER_DOMAIN) { /* Terminate corner drag mode */ glDisable(GL_COLOR_LOGIC_OP); } mode = SWAP_EDGES; glutMouseFunc(swapEdges); glutSetCursor(GLUT_CURSOR_CROSSHAIR); /* Change cursor */ return; case 'f': /* Fill R with triangles */ fillRt(); break; case 'F': /* Write the triangle mesh to the output file */ trwrite(fname, nc, nb, nv, nt, ne, vtxy, uv, ltri); break; case 'H': /* Toggle use of off-center/circumcenter in refineRd */ offcenter = !offcenter; return; case 'i': /* Triangulation integrity check */ trmtst(nc, nb, nv, nn, nt, ne, vtxy, ltri, lvtx, ledg, lnrv, lcrv, listn); return; case 'j': /* Toggle antialiasing */ antialias = !antialias; if (antialias) { glHint(GL_LINE_SMOOTH_HINT, GL_NICEST); glEnable(GL_LINE_SMOOTH); glEnable(GL_BLEND); glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); } else { glDisable(GL_LINE_SMOOTH); glDisable(GL_BLEND); } break; case 'k': /* Decrease line width */ line_width *= 0.80; break; case 'K': /* Increase line width */ line_width *= 1.25; break; case 'l': /* Get a new title location title_pos */ mode = PLACE_TITLE; glutMouseFunc(getTitlePos); glutSetCursor(GLUT_CURSOR_RIGHT_ARROW); /* Change cursor */ break; case 'L': /* Toggle boundary layer smoothing */ bdry_layer = !bdry_layer; return; case 'm': /* Decrease minimum angle (1 to 35) for refineRd */ if (min_angle > 1) min_angle--; c = cos(acos(-1.0)*min_angle/180.0); scamax = c*c; offscale = 0.48*sqrt((1.0+c)/(1.0-c)); break; case 'M': /* Increase minimum angle (in degrees) for refineRd */ if (min_angle < 35) min_angle++; c = cos(acos(-1.0)*min_angle/180.0); scamax = c*c; offscale = 0.48*sqrt((1.0+c)/(1.0-c)); break; case 'n': /* Convert vertices in R to removable */ convertRv(1); break; case 'N': /* Convert vertices in R to non-removable */ convertRv(0); break; case 'o': /* Draw a circle by dragging a point on the circumference with the left mouse button */ glutMouseFunc(getCircle); glutSetCursor(GLUT_CURSOR_CROSSHAIR); /* Change cursor */ break; case 'O': /* Project non-removable vertices in R onto circle */ moveRn(); break; case 'p': /* Pan: select center (xc,yc) */ mode = PAN; glutMouseFunc(getCenter); glutSetCursor(GLUT_CURSOR_CROSSHAIR); /* Change cursor */ break; case 'P': /* Print triangle mesh data structure */ trprint(nc, nb, nv, nn, nt, ne, vtxy, uv, ltri, lvtx, ledg, lnrv, lcrv); return; case 'q': /* Query current state */ printf("\n\n"); printf(" Mesh parameters: nc = %u, nb = %u, nv = %u, nn = %u" ", nt = %u, ne = %u\n", nc, nb, nv, nn, nt, ne); printf(" nvmax = %u\n\n", nvmax); printf(" View volume parameters:\n" " Center x = %g\n", xc); printf(" Center y = %g\n", yc); printf(" Width = %g\n", dx); printf(" Height = %g\n\n", dy); printf(" Boundary layer flag = %s\n", flagName[bdry_layer]); printf(" Antialiasing flag = %s\n", flagName[antialias]); printf(" Preserve aspect ratio flag = %s\n", flagName[aspect]); printf(" Offcenter flag = %s\n", flagName[offcenter]); printf(" Line width = %.3f\n", line_width); printf(" Font number for indices = %u\n", ind_font); printf(" Font number for title = %u\n", title_font); printf(" Minimum angle for refineRd = %d\n", min_angle); printf(" Maximum number of Steiner points per refineRd = %d\n", nspmax); printf(" Number of vertices in boundary of R = %d\n", np); printf(" Current mode = %s\n", modeName[mode]); return; case 'r': /* Restore default viewing volume */ dx = 1.2*(xmax-xmin); dy = 1.2*(ymax-ymin); xc = (xmin+xmax)/2.0; yc = (ymin+ymax)/2.0; reshape(glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT)); break; case 'R': /* Select a new polygon R */ mode = CONSTRUCT_R; newR = 1; np = 0; glutMouseFunc(getPolygon); glutSetCursor(GLUT_CURSOR_CROSSHAIR); /* Change cursor */ break; case 's': /* Optimize R with edge swaps */ optimizeRt(); break; case 'S': /* Optimize R by applying edge swaps, and then smoothing (moving vertices) */ optimizeRt(); optimizeRv(); break; case 't': /* Toggle display of triangle indices */ plot_t = !plot_t; break; case 'T': /* Input a new character string for title */ printf("\n Enter a title string: "); getLin(title, LSTRING); break; case 'u': /* Enter uv-edit mode */ if (mode == CONSTRUCT_R) { /* Terminate polygon selection */ glDisable(GL_COLOR_LOGIC_OP); np = 0; } if (mode == ALTER_DOMAIN) { /* Terminate corner drag mode */ glDisable(GL_COLOR_LOGIC_OP); } mode = EDIT_UV; glutMouseFunc(getVertex); glutSetCursor(GLUT_CURSOR_CROSSHAIR); /* Change cursor */ return; case 'v': /* Toggle display of vertex indices */ plot_v = !plot_v; break; case 'V': /* Toggle display of non-removable vertex indices */ plot_nrv = !plot_nrv; break; case 'w': /* Select viewing volume center, width, and height */ glutMouseFunc(getViewVolume); glutSetCursor(GLUT_CURSOR_CROSSHAIR); /* Change cursor */ return; case 'x': /* Halve nspmax */ nspmax *= 0.5; return; case 'X': /* Double nspmax */ nspmax *= 2.0; return; case 'z': /* Zoom in */ dx *= 0.8; dy *= 0.8; reshape(glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT)); break; case 'Z': /* Zoom out */ dx *= 1.25; dy *= 1.25; reshape(glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT)); break; case '&': /* Refine R using Delaunay refinement */ refineRd(); break; case '+': /* Refine R by adding triangle barycenters */ refineRt(); break; case '-': /* Unrefine the mesh by reversing the previous refinement by refineRe or refineRt */ unrefine(); break; case '>': /* Refine R by adding edge midpoints */ refineRe(); break; case '<': /* Coarsen R by removing interior vertices */ coarsenRv(); break; case 8: /* Remove triangles in R */ case 127: /* Remove triangles in R */ removeRt(); break; case '!': /* Toggle font selection ('0' to '8') between indices (font_select = 0) and title (font_select = 1) */ font_select = !font_select; return; case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': if (!font_select) ind_font = item - '0'; else title_font = item - '0'; break; case 10: /* Terminate current mode */ case 13: if (mode == CONSTRUCT_R) closeR(); mode = DEFAULT; /* Restore default mode */ glutMouseFunc(mouse); /* Restore default mouse callback */ glutSetCursor(GLUT_CURSOR_INHERIT); /* Standard cursor */ glDisable(GL_COLOR_LOGIC_OP); /* Overwrite mode */ break; case 27: /* Terminate program or current mode */ if (mode == CONSTRUCT_R) { mode = DEFAULT; /* Terminate polygon selection mode */ glutMouseFunc(mouse); glutSetCursor(GLUT_CURSOR_INHERIT); glDisable(GL_COLOR_LOGIC_OP); np = 0; break; } if (mode == ALTER_DOMAIN || mode == SWAP_EDGES || mode == EDIT_UV) { mode = DEFAULT; /* Terminate drag mode */ glutMouseFunc(mouse); glutSetCursor(GLUT_CURSOR_INHERIT); glDisable(GL_COLOR_LOGIC_OP); break; } sprintf(dlg_string,"You really want to quit?"); dlg_exit = GL_TRUE; glutHideWindow(); /* Hide window 0 */ glutSetWindow(window[1]); /* Make window 1 current */ glutShowWindow(); /* Make window 1 visible */ dlgDisplay(); break; default: break; } glutPostRedisplay(); return; } /* End of menu */ /* *-------------------------------------------------------------------- * * motion: Callback triggered by mouse motion with one or more mouse * buttons pressed. * * If mode = CONSTRUCT_R, this function, operating in XOR write mode, * draws the line segment with endpoints (xv1,yv1) and (xv2,yv2) * (thus erasing the previously drawn line segments), updates (xv2, * yv2) to the object coordinates to the mouse position, and draws * the updated line segments. * * If mode = ALTER_R, this function, operating in XOR write mode, * draws the two line segments with endpoints (xv1,yv1), (xv2,yv2), * and (xv3,yv3) (thus erasing the previously drawn line segments), * updates (xv2,yv2) to the object coordinates to the mouse position, * and draws the updated line segments. * * If mode = ALTER_VIEW_VOLUME, this function, operating in XOR write * mode, draws an outline of the rectangle [xv1,xv2] X [yv1,yv2] * (thus erasing the previously drawn rectangle), updates (xv2,yv2) * to the object coordinates corresponding to the mouse position, and * draws an outline of the updated rectangle. * * If mode = SELECT_CIRCLE, this function, operating in XOR write * mode, draws a circle of radius cr centered at (cx,cy), thus * erasing the previously drawn circle, changes cr to the distance * from (cx,cy) to the mouse position, and draws the circle with the * new radius. * * If nnbv > 1 (and mode = ALTER_DOMAIN), this function, operating in * XOR write mode, draws the set of line segments from (xv1,yv1) to * the neighbors of vertex kv1, erasing the previously drawn line * segments, updates (xv1,yv1) to the mouse position, and draws the * new set of line segments. The indices of the neighbors of kv1 * are stored in the first nnbv positions of lste. * * glmesh2 functions called: displayPosition, drawCircle, xlate * *-------------------------------------------------------------------- */ static void motion(int x, int y) { GLdouble xm, ym, zm; /* Computed object coordinates */ GLdouble x2, y2; GLint i, kv, nnb; xlate(x,y,0.0, &xm,&ym,&zm); if (mode == CONSTRUCT_R) { /* * Erase line segment from (xv1,yv1) to (xv2,yv2). */ glBegin(GL_LINES); glVertex2d(xv1, yv1); glVertex2d(xv2, yv2); glEnd(); /* * Update (xv2,yv2). */ xv2 = xm; yv2 = ym; /* * Draw line segment from (xv1,yv1) to (xv2,yv2). */ glBegin(GL_LINES); glVertex2d(xv1, yv1); glVertex2d(xv2, yv2); glEnd(); } else if (mode == ALTER_R) { /* * Erase line segments from (xv1,yv1) to (xv2,yv2) and (xv2,yv2) to * (xv3,yv3). */ glBegin(GL_LINE_STRIP); glVertex2d(xv1, yv1); glVertex2d(xv2, yv2); glVertex2d(xv3, yv3); glEnd(); /* * Update (xv2,yv2). */ xv2 = xm; yv2 = ym; /* * Draw line segments from (xv1,yv1) to (xv2,yv2) and (xv2,yv2) to * (xv3,yv3). */ glBegin(GL_LINE_STRIP); glVertex2d(xv1, yv1); glVertex2d(xv2, yv2); glVertex2d(xv3, yv3); glEnd(); } else if (mode == ALTER_VIEW_VOLUME) { /* * Erase rectangle [xv1,xv2] X [yv1,yv2]. */ glBegin(GL_LINE_LOOP); glVertex2d(xv1, yv1); glVertex2d(xv2, yv1); glVertex2d(xv2, yv2); glVertex2d(xv1, yv2); glEnd(); /* * Update (xv2,yv2). */ xv2 = xm; yv2 = ym; /* * Draw rectangle [xv1,xv2] X [yv1,yv2]. */ glBegin(GL_LINE_LOOP); glVertex2d(xv1, yv1); glVertex2d(xv2, yv1); glVertex2d(xv2, yv2); glVertex2d(xv1, yv2); glEnd(); } else if (mode == SELECT_CIRCLE) { /* * mode = SELECT_CIRCLE. */ drawCircle(); cr = sqrt( (xm-cx)*(xm-cx) + (ym-cy)*(ym-cy) ); drawCircle(); } else if (nnbv > 1) { /* * mode = ALTER_DOMAIN and nnbv > 1. */ nnb = nnbv; glBegin(GL_LINES); for (i = 0; i < nnb; i++) { kv = lste[i]; x2 = vtxy[kv][0]; y2 = vtxy[kv][1]; glVertex2d(xv1, yv1); glVertex2d(x2, y2); } glEnd(); xv1 = xm; yv1 = ym; glBegin(GL_LINES); for (i = 0; i < nnb; i++) { kv = lste[i]; x2 = vtxy[kv][0]; y2 = vtxy[kv][1]; glVertex2d(xv1, yv1); glVertex2d(x2, y2); } glEnd(); } /* * Display current mouse position in the upper left corner. */ displayPosition(xm,ym); return; } /* End of motion */ /* *-------------------------------------------------------------------- * * mouse: Callback triggered by a mouse button press or release with * the mouse position in the primary window (except for * buttons attached to menus). * * A left button press may be used to select and drag a vertex of the * polygon boundary R, insert a new vertex in R by selecting a * boundary point other than a vertex, or delete a vertex by * selecting it with the shift key engaged. * * glmesh2 functions called: intsec, plDistance, xlate * *-------------------------------------------------------------------- */ static void mouse(int button, int state, int x, int y) { GLdouble wx, wy, wz; /* Computed object coordinates */ GLdouble tols; GLint i, ip1, k, km1, km2, kp1; GLboolean intr; if (button != GLUT_LEFT_BUTTON || np == 0) return; if (state == GLUT_DOWN) { /* * Compute the object coordinates (wx,wy) of the mouse cursor. */ xlate(x,y,0.0, &wx,&wy,&wz); /* * Find a vertex R[k], if any, within squared distance tols of the * mouse position (wx,wy). */ tols = (sx <= sy) ? sx : sy; /* tols = min(sx, sy) */ tols = 25.0*tols*tols; /* 5-pixel tolerance */ intr = 0; for (k = 0; k <= np-1; k++) { if ((intr = ((wx-R[k][0])*(wx-R[k][0]) + (wy-R[k][1])*(wy-R[k][1]) <= tols))) break; } if (glutGetModifiers() == GLUT_ACTIVE_SHIFT) { if (!intr || np < 4) { return; } else { /* * Vertex selected with shift key down, and np > 3: delete vertex k. */ for (i = k; i < np-1; i++) { R[i][0] = R[i+1][0]; R[i][1] = R[i+1][1]; } np--; newR = 1; glutPostRedisplay(); return; } } if (!intr) { /* * No shift key and no vertex selected. Test for the mouse position * within a tolerance of the polygon boundary. */ tols = 0.64*tols; /* 4-pixel tolerance */ intr = 0; for (k = 0; k < np; k++) { kp1 = (k+1)%np; if ((intr = plDistance(wx, wy, R[k][0], R[k][1], R[kp1][0], R[kp1][1], tols))) break; } if (intr) { /* * Insert a new vertex after reallocating memory for R if necessary. */ if (np >= npmax) { npmax *= 2; R = realloc(R, npmax*sizeof *R); if (R == NULL) { printf("glmesh2: Unable to allocate sufficient " "memory for R.\n"); exit(1); } } for (i = np-1; i > k; i--) { R[i+1][0] = R[i][0]; R[i+1][1] = R[i][1]; } k++; R[k][0] = wx; R[k][1] = wy; np++; } else { return; } } /* * Set up for dragging (Function motion): store R[k-1], R[k], and * R[k+1] in (xv1,yv1), (xv2,yv2), and (xv3,yv3), set XOR write mode, * set newR = 1, mode = ALTER_R, set indxR = k, and disable the cursor. */ km1 = (k-1+np)%np; xv1 = R[km1][0]; yv1 = R[km1][1]; xv2 = R[k][0]; yv2 = R[k][1]; kp1 = (k+1)%np; xv3 = R[kp1][0]; yv3 = R[kp1][1]; glEnable(GL_COLOR_LOGIC_OP); glLogicOp(GL_XOR); newR = 1; mode = ALTER_R; indxR = k; glutSetCursor(GLUT_CURSOR_NONE); } else { /* * Left button release: test for an intersection. */ if (mode != ALTER_R) return; k = indxR; km1 = (k-1+np)%np; km2 = (k-2+np)%np; intr = 0; for (i = 0; i < np; i++) { ip1 = (i+1)%np; if (i == km2 || i == km1 || i == k) continue; if (intsec(xv1,yv1,xv2,yv2,R[i][0],R[i][1],R[ip1][0], R[ip1][1])) { intr = 1; break; } } if (!intr && np > 3) { intr = intsec(xv2,yv2,xv3,yv3,R[km2][0],R[km2][1],R[km1][0], R[km1][1]); } if (intr) { /* * Invalid choice: erase the line segments. */ glBegin(GL_LINE_STRIP); glVertex2d(xv1,yv1); glVertex2d(xv2,yv2); glVertex2d(xv3,yv3); glEnd(); glutSwapBuffers(); } else { /* * Update R. */ R[k][0] = xv2; R[k][1] = yv2; } /* * Restore the inherited cursor, and disable XOR write mode. */ mode = DEFAULT; glutSetCursor(GLUT_CURSOR_INHERIT); glDisable(GL_COLOR_LOGIC_OP); glutPostRedisplay(); } return; } /* End of mouse */ /* *-------------------------------------------------------------------- * * moveRn: Project all non-removable vertices in R (that can be * moved without destroying the mesh) onto the current * circle (if any). * * glmesh2 functions called: bbox, initR, intsec * *-------------------------------------------------------------------- */ static void moveRn(void) { double s, x1, x2, x3, xp, y1, y2, y3, yp; int i, k, kt, kv, kv1, kv2, kv3 = 0, kvn; unsigned char intr, visible; /* * Test for a circle. */ if (cr <= 0.0) return; /* * Store listR if necessary. */ if (newR) initR(); /* * Outer loop on non-removable vertices kv1 = lnrv[k] in R. */ for (k = 0; k < nn; k++) { kv1 = lnrv[k]; if (!listR[kv1]) continue; /* * Compute the coordinates (xp,yp) of the projection p of kv1 * onto the circle. If kv1 is not at the center (cx,cy), it can * be projeced onto the circle (but it could end up anywhere on * the circle if it starts too near the center). */ x1 = vtxy[kv1][0]; y1 = vtxy[kv1][1]; s = sqrt( (x1-cx)*(x1-cx) + (y1-cy)*(y1-cy) ); if (s <= 0.0) continue; s = cr/s; xp = cx + s*(x1-cx); yp = cy + s*(y1-cy); /* * The new point p is visible iff p Left kv2 -> kv3 for all pairs of * adjacent neighbors (kv2,kv3) of kv1, and there is no intersection * of the path p-kv1 with a boundary edge not containing kv1 in the * same boundary curve. */ visible = 1; kt = lvtx[kv1]; while (kt >= 0) { if (ltri[kt][0] == kv1) i = 1; else if (ltri[kt][1] == kv1) i = 2; else i = 0; kv2 = ltri[kt][i]; kv3 = ltri[kt][(i+1)%3]; x2 = vtxy[kv2][0]; y2 = vtxy[kv2][1]; x3 = vtxy[kv3][0]; y3 = vtxy[kv3][1]; if ( (x3-x2)*(yp-y2) < (y3-y2)*(xp-x2) ) { visible = 0; break; } kt = ltri[kt][i+3]; } if (visible) { /* * Set kv2 and kv3 to the first and last neighbors of kv1, and * test for intersections. */ intr = 0; kt = lvtx[kv1]; if (ltri[kt][0] == kv1) i = 1; else if (ltri[kt][1] == kv1) i = 2; else i = 0; kv2 = ltri[kt][i]; kv = kv2; while (kv != kv3) { kt = lvtx[kv]; if (ltri[kt][0] == kv) i = 1; else if (ltri[kt][1] == kv) i = 2; else i = 0; kvn = ltri[kt][i]; if (intsec(xp,yp,x1,y1,vtxy[kv][0],vtxy[kv][1], vtxy[kvn][0],vtxy[kvn][1])) { intr = 1; break; } kv = kvn; } if (intr) visible = 0; } if (visible) { /* * Move kv1 to its new position. */ vtxy[kv1][0] = xp; vtxy[kv1][1] = yp; } } /* * Update the bounding box. */ bbox(); return; } /* End of moveRn */ /* *-------------------------------------------------------------------- * * offCenter: Compute an off-center of a triangle: a point on the * perpendicular bisector of the shortest edge, between * the midpoint of the edge and the circumcenter of the * triangle with distance from the edge just small enough * that the new triangle formed by the edge and the off- * center will not be split. The off-center is the * circumcenter if the circumcenter is close enough to * the edge or if global variable offcenter = FALSE (and * the the circumcenter is computable). * * On input: kt is the index (0 to nt-1) of the triangle whose * off-center is to be computed, and imin is the local * index (0, 1, or 2) of the vertex opposite the shortest * edge of triangle kt. * * On output: ctrx and ctry are the Cartesian coordinates of the * off-center (or circumcenter). * * Reference: Alper Ungor, Off-centers: A new type of Steiner * points for computing size-optimal quality-guaranteed * Delaunay triangulations, Computational Geometry: * Theory and Applications, 42 (2), Feb 2009, 109-118. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void offCenter(int kt, int imin, double *ctrx, double *ctry) { double ds1, ds2, dx1, dx2, dxc, dxo, dy1, dy2, dyc, dyo, s, x0, x1, x2, y0, y1, y2; int kv0, kv1, kv2; /* * Compute the coordinates of the vertices v0, v1, v2, where v0 and * v1 are the endpoints of the shortest edge. All vectors are * computed relative to origin v0. */ kv2 = ltri[kt][imin]; kv0 = ltri[kt][(imin+1)%3]; kv1 = ltri[kt][(imin+2)%3]; x0 = vtxy[kv0][0]; y0 = vtxy[kv0][1]; x1 = vtxy[kv1][0]; y1 = vtxy[kv1][1]; x2 = vtxy[kv2][0]; y2 = vtxy[kv2][1]; dx1 = x1 - x0; dy1 = y1 - y0; dx2 = x2 - x0; dy2 = y2 - y0; /* * s = (v1-v0) X (v2-v0) */ s = dx1*dy2 - dx2*dy1; if (s <= 0.0) { printf("\n offCenter: triangle %d has non-positive signed " "area.\n", kt); *ctrx = x0 + 0.5*dx1 - offscale*dy1; *ctry = y0 + 0.5*dy1 + offscale*dx1; return; } s = 0.5/s; /* * ds1 = |v1-v0|^2, ds2 = |v2-v0|^2, and * (dxc,dyc) = c - v0, where c is the circumcenter. */ ds1 = dx1*dx1 + dy1*dy1; ds2 = dx2*dx2 + dy2*dy2; dxc = s*(dy2*ds1 - dy1*ds2); dyc = s*(dx1*ds2 - dx2*ds1); /* * Overwrite (dxc,dyc) with o - v0, for off-center o, if * offcenter = TRUE, offscale > 0 (min_angle > 0), and the * off-center is between the edge and the circumcenter. */ if (offcenter && offscale > 0.0) { dxo = 0.5*dx1 - offscale*dy1; dyo = 0.5*dy1 + offscale*dx1; if (dxo*dxo + dyo*dyo < dxc*dxc + dyc*dyc) { dxc = dxo; dyc = dyo; } } *ctrx = x0 + dxc; *ctry = y0 + dyc; return; } /* End of offCenter */ /* *-------------------------------------------------------------------- * * optimizeRt: Optimize R with Delaunay edge swaps in convex * quadrilaterals consisting of pairs of adjacent * triangles in R. * * The indices of the interior edges associated with pairs of * adjacent triangles contained in R are stored in lste. Then a * sequence of sweeps is executed in which each sweep consists of * applying the swap test (empty circumcircle test) and appropriate * swaps to the interior edges in the order in which they are stored. * The iteration is repeated until no swap occurs or maxit iterations * have been performed. The bound on the number of iterations * (defined below) may be necessary to prevent an infinite loop * caused by cycling (reversing the effect of a previous swap) due to * floating point inaccuracy when four or more vertices are nearly * cocircular. * * glmesh2 functions called: initR, swap, swptst * *-------------------------------------------------------------------- */ static void optimizeRt(void) { int maxit = 4; /* Maximum number of iterations */ int swp; /* Count of number of swaps in a pass */ int in1, in2, io1, io2; int i, iter, j, ke, kt, ktn, le; /* * Store listR if necessary. */ if (newR) initR(); /* * Store the interior edge indices in lste with count le. * Each pair of adjacent triangles in R is encountered twice * (once for each triangle). The corresponding edge is only * stored once by associating it with the smaller triangle index. */ le = 0; for (kt = 0; kt < nt; kt++) { if (!listR[ltri[kt][0]] || !listR[ltri[kt][1]] || !listR[ltri[kt][2]]) continue; for (i = 0; i < 3; i++) { ktn = ltri[kt][i+3]; if (ktn < kt) continue; if (listR[ltri[ktn][0]] && listR[ltri[ktn][1]] && listR[ltri[ktn][2]]) { lste[le] = ltri[kt][i+6]; le++; } } } if (le == 0) return; /* * Top of outer loop. */ for (iter = 0; iter < maxit; iter++) { swp = 0; /* * Inner loop on interior edges. */ for (j = 0; j < le; j++) { ke = lste[j]; kt = ledg[ke]; if (ltri[kt][6] == ke) i = 0; else if (ltri[kt][7] == ke) i = 1; else i = 2; in1 = ltri[kt][i]; io1 = ltri[kt][(i+1)%3]; io2 = ltri[kt][(i+2)%3]; kt = ltri[kt][i+3]; if (ltri[kt][6] == ke) i = 0; else if (ltri[kt][7] == ke) i = 1; else i = 2; in2 = ltri[kt][i]; if (swptst(in1,in2,io1,io2)) { swap(ke); swp++; } } if (!swp) break; } return; } /* End of optimizeRt */ /* *-------------------------------------------------------------------- * * optimizeRv: Optimize R by applying a single iteration of a CVT * (Centroid Voronoi Tessellation) smoothing operation: * sweep through the removable vertices, moving each * boundary vertex to the midpoint of its predecessor * and successor (unless they are adjacent), and each * interior vertex to a weighted sum (convex combina- * tion) of the centroids of the triangles that contain * the vertex, where the weights are the triangle areas * relative to the star area. This is similar to * Laplacian smoothing by a Gauss-Seidel iteration. If * global variable bdry_layer = TRUE, each interior * vertex adjacent to a boundary vertex is moved to a * location that depends only on the triangles that * contain the vertex and contain a boundary vertex, * so that the new location is usually closer to the * boundary. * * Reference: Long Chen, Mesh smoothing schemes based on optimal * Delaunay triangulation, in Proceedings of 13th * International Meshing Roundtable, 2004, 109-120. * * glmesh2 functions called: initR, intsec * *-------------------------------------------------------------------- */ static void optimizeRv(void) { double s, sumw, tolc, tols, w, x, x0, x1, x2, y, y0, y1, y2; unsigned char bdry, bdry1, bdry2, visible; int i, j, kt, kt1, kv0, kv1, kv2, kvn, nnb; /* * Store a tolerance used by the left test to ensure visibility of a * new interior vertex position. */ tolc = 1.e-4; /* * Compute a tolerance tols on the squared length of a boundary * edge chosen to make the relative length greater than 1.e-3. */ tols = 1.e-3*dx; tols = tols*tols; /* * Store listR if necessary. */ if (newR) initR(); /* * Loop on vertices kv0 in R. */ for (kv0 = 0; kv0 < nv; kv0++) { if (!listR[kv0]) continue; x0 = vtxy[kv0][0]; y0 = vtxy[kv0][1]; /* * Test for kv0 removable. */ if (listn[kv0]) continue; /* * Store the neighbors of kv0 in the first nnb positions of lste, * where the first neighbor is also the last neighbor (stored twice) * unless kv0 is a boundary vertex. */ nnb = 0; kt1 = lvtx[kv0]; kt = kt1; do { if (ltri[kt][0] == kv0) i = 1; else if (ltri[kt][1] == kv0) i = 2; else i = 0; lste[nnb] = ltri[kt][i]; nnb++; kvn = ltri[kt][(i+1)%3]; kt = ltri[kt][i+3]; } while (kt >= 0 && kt != kt1); lste[nnb] = kvn; nnb++; visible = 1; if (kt < 0) { /* * Boundary vertex: compute the new position (x,y) as the midpoint * of the first and last neighbors kv1 and kv2. */ kv1 = lste[0]; kv2 = lste[nnb-1]; x1 = vtxy[kv1][0]; y1 = vtxy[kv1][1]; x2 = vtxy[kv2][0]; y2 = vtxy[kv2][1]; x = 0.5*(x1 + x2); y = 0.5*(y1 + y2); /* * The boundary is not altered if the squared length of edge kv1-kv2 * is less than tols. This test prevents a hole defined by an * interior boundary curve from shrinking to a point, causing * numerical problems. */ if ( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) < tols ) continue; } else { /* * Interior vertex: set bdry = TRUE iff bdry_layer = TRUE and a * neighbor kv1 of kv0 is a boundary vertex. */ bdry = 0; if (bdry_layer) { for (i = 0; i < nnb-1; i++) { kv1 = lste[i]; kt = lvtx[kv1]; if (ltri[kt][0] == kv1) j = 2; else if (ltri[kt][1] == kv1) j = 0; else j = 1; if ((bdry = (ltri[kt][j+3] < 0))) break; } } /* * Compute the potential new position as a weigted sum (x,y) of the * centroids of the triangles (kv0,kv1,kv2). * * Initialization: */ x = 0.0; y = 0.0; sumw = 0.0; bdry1 = 0; bdry2 = 0; if (bdry) { kv1 = lste[0]; kt = lvtx[kv1]; if (ltri[kt][0] == kv1) j = 2; else if (ltri[kt][1] == kv1) j = 0; else j = 1; bdry2 = (ltri[kt][j+3] < 0); } for (i = 0; i < nnb-1; i++) { bdry1 = bdry2; kv1 = lste[i]; kv2 = lste[i+1]; if (bdry) { /* * Bypass the triangle if neither kv1 nor kv2 is a boundary vertex. */ kt = lvtx[kv2]; if (ltri[kt][0] == kv2) j = 2; else if (ltri[kt][1] == kv2) j = 0; else j = 1; bdry2 = (ltri[kt][j+3] < 0); if (!bdry1 && !bdry2) continue; } x1 = vtxy[kv1][0]; y1 = vtxy[kv1][1]; x2 = vtxy[kv2][0]; y2 = vtxy[kv2][1]; w = (x2-x1)*(y0-y1) - (y2-y1)*(x0-x1); sumw += w; x += w*(x0 + x1 + x2); y += w*(y0 + y1 + y2); } /* * Abort if sumw = 0 (all vertices in the star are collinear). */ if (!sumw) continue; /* * Scale (x,y). */ s = 1.0/(3.0*sumw); x *= s; y *= s; } /* * Test for (x,y) visible from the neighbors: (x,y) Strictly Left * kv1 -> kv2 for all pairs of adjacent neighbors (kv1,kv2). A * boundary vertex with two neighbors is a special case: (x,y) is * the midpoint of kv1 -> kv2, and is necessarily visible. */ if (nnb > 2) { for (i = 0; i < nnb-1; i++) { kv1 = lste[i]; kv2 = lste[i+1]; x1 = vtxy[kv1][0]; y1 = vtxy[kv1][1]; x2 = vtxy[kv2][0]; y2 = vtxy[kv2][1]; if ( (x2-x1)*(y-y1) - (y2-y1)*(x-x1) < tolc ) { visible = 0; break; } } } if (visible) { /* * Move kv0 to its new position. */ vtxy[kv0][0] = x; vtxy[kv0][1] = y; } } return; } /* End of optimizeRv */ /* *-------------------------------------------------------------------- * * passiveMotion: Callback triggered by mouse motion with no mouse * buttons pressed. * * This function displays the object coordinates of the mouse * position. * * glmesh2 functions called: displayPosition, xlate * *-------------------------------------------------------------------- */ static void passiveMotion(int x, int y) { GLdouble xm, ym, zm; /* Computed object coordinates */ /* * Compute the object coordinates (xm,ym) of the mouse cursor. */ xlate(x,y,0.0, &xm,&ym,&zm); /* * Display current mouse position in the upper left corner. */ displayPosition(xm,ym); return; } /* End of passiveMotion */ /* *-------------------------------------------------------------------- * * plDistance: Compute a point-line distance. The return value is * TRUE iff the projection p* of point p onto the line * defined by p1 and p2 lies on the line segment p1-p2 * and satisfied <= tols. * * On input: p = (xp,yp), p1 = (x1,y1), p2 = (x2,y2), and tols is * the maximum squared distance between p and p* for which * the return value is TRUE. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static unsigned char plDistance(double xp, double yp, double x1, double y1, double x2, double y2, double tols) { double cp, dp, ds12, dx12, dx1p, dy12, dy1p; /* * p* = p1 + s*(p2-p1), where = 0, so that * s = dp/ds12 for dp = and ds12 = . * Then p* lies on the line segment iff 0 <= s <= 1 for ds12 > 0. * The return value is FALSE if p1 = p2. */ dx12 = x2-x1; dy12 = y2-y1; ds12 = dx12*dx12 + dy12*dy12; dx1p = xp-x1; dy1p = yp-y1; dp = dx1p*dx12 + dy1p*dy12; if (dp < 0.0 || dp > ds12) return 0; /* * The squared distance is ds = = cp*cp/ds12 for * cross product (z-component) cp = (p2-p1) X (p-p1). */ cp = dx12*dy1p - dy12*dx1p; return (cp*cp <= tols*ds12); } /* End of plDistance */ /* *-------------------------------------------------------------------- * * reallocQe: Reallocate storage for Qe, and copy the contents of * the old array into the new array with an ordering * that reflects the circular array storage scheme. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void reallocQe(void) { int i, ip, n0; struct boundary_edge *ptr; /* * Save the old queue size and pointer in n0 and ptr, double neqmax, * and allocate a new array. */ n0 = neqmax; ptr = Qe; neqmax *= 2; Qe = malloc(neqmax*sizeof *Qe); /* * Copy the contents of the queue to the new addresses. */ ip = qefront; for (i = 0; i < n0; i++) { Qe[i].index_e = ptr[ip].index_e; Qe[i].index_v1 = ptr[ip].index_v1; Qe[i].index_v2 = ptr[ip].index_v2; ip = (ip+1)%n0; } /* * Update the front and back pointers, and free the old memory. */ qefront = 0; qeback = n0-1; free(ptr); return; } /* End of reallocQe */ /* *-------------------------------------------------------------------- * * refineRd: Refine R by Delaunay refinement. Bad triangles are * split by inserting their circumcenters, and encroached * boundary edges are split by inserting their midpoints. * Edges are then swapped as required by the empty circum- * circle criterion (Function swptst). If a constrained * Delaunay triangulation is input, a Delaunay triangula- * tion will be output. * * The circumcenter of a triangle may not be contained in the * triangle, but its insertion into a constrained Delaunay triangu- * lation will cause at least one of its edges to be swapped out. * A boundary edge is said to be encroached upon by a vertex in its * diametral circle (the circle that has the edge as a diameter). * If a circumcenter encroaches upon a boundary edge, it is not * inserted. Note that, if a circumcenter is not contained in a * triangle, it must encroach upon a boundary edge, and is therefore * not inserted. * * An encroached boundary edge may be split at a point other than the * exact midpoint. Refer to Function splitEdges. * * If offcenter = TRUE, the off-center is used in place of the * circumcenter (Function offCenter). * * References: * * Jim Rupert, A Delaunay refinement algorithm for quality 2- * dimensional mesh generation, Journal of Algorithms 18 (3), * May 1995, 548-585. * * Jonathan Richard Shewchuk, Delaunay refinement algorithms for * triangular mesh generation, Computational Geometry: Theory * and Applications 22 (1-3), May 2002, 21-74. * * glmesh2 functions called: buildQt, deleteQt, delvtx, encroached, * initR, insertQt, insertv, offCenter, * optimizeRt, splitEdges, trealloc, * trfind, unswap, update * *-------------------------------------------------------------------- */ static void refineRd(void) { double b1, b2, b3, ctrx, ctry, emin; int i, imin, k, ke, kt, kt0, kv1, kv2, kv3, neq0 = 0, ns; struct triangle ptri; /* * Store listR if necessary. */ if (newR) initR(); /* * Initialize nsp, and push nv onto rstack. */ nsp = 0; if (nref == nrefmax) { nrefmax *= 2; rstack = realloc(rstack, nrefmax*sizeof *rstack); if (rstack == NULL) { printf("glmesh2: Unable to allocate sufficient memory " "for rstack.\n"); exit(1); } } rstack[nref] = nv; nref++; /* * Apply the necessary edge swaps to create a constrained Delaunay * triangulation. */ optimizeRt(); /* * Allocate storage for a queue Qe of boundary edges implemented as * a circular array: front and back indices, qefront and qeback, * are incremented modulo neqmax for dequeueing and enqueueing, * respectively. Note that qeback = (qefront-1)%neqmax in both the * case of an empty queue (neq = 0) and a full queue (neq = neqmax). */ neqmax = nb; Qe = malloc(neqmax*sizeof *Qe); /* * Enqueue all encroached boundary edges in R, and then split them. */ neq = 0; qefront = 0; qeback = -1; for (ke = 0; ke < ne; ke++) { kt = ledg[ke]; if (ltri[kt][6] == ke) i = 0; else if (ltri[kt][7] == ke) i = 1; else i = 2; if (ltri[kt][i+3] >= 0) continue; kv1 = ltri[kt][(i+1)%3]; kv2 = ltri[kt][(i+2)%3]; if (!listR[kv1] || !listR[kv2]) continue; if (encroached(ke, &kv1, &kv2)) { if (neq == neqmax) reallocQe(); qeback = (qeback+1)%neqmax; Qe[qeback].index_e = ke; Qe[qeback].index_v1 = kv1; Qe[qeback].index_v2 = kv2; neq++; } } splitEdges(0); /* * Test for failure to converge. */ if (neq) printf("\n splitEdges aborted with neq = %d\n",neq); /* * Construct the priority queue of bad triangles. */ buildQt(); /* * Loop on triangles in Qt. */ while (ntq && nsp <= nspmax) { /* * Dequeue the next triangle kt0. */ deleteQt(&ptri); kt0 = ptri.index_t; kv1 = ptri.index_v1; kv2 = ptri.index_v2; kv3 = ptri.index_v3; imin = ptri.imin; emin = ptri.emin; /* * Bypass triangles that have been removed from the triangulation. */ if (ltri[kt0][0] != kv1 || ltri[kt0][1] != kv2 || ltri[kt0][2] != kv3) continue; /* * Compute the off-center (ctrx,ctry) of triangle kt0. */ offCenter(kt0, imin, &ctrx, &ctry); /* * Find the triangle kt, if any, that contains (ctrx,ctry), along * with the barycentric coordinates b1, b2, and b3. If a boundary * edge intersects the path from the interior of kt0 to (ctrx,ctry), * its index is returned in ke. The return value is 2 if (ctrx,ctry) * is found to lie in a null triangle. */ ke = -1; i = trfind(kt0, ctrx, ctry, &ke, &kt, &b1, &b2, &b3); if (i) { if (b1 < 0.0 || b2 < 0.0 || b3 < 0.0 || b1 != b1 || b2 != b2 || b3 != b3) { printf("\n trfind: b1 = %e, b2 = %e, b3 = %e\n", b1, b2, b3); break; } if (i > 1) break; /* * No boundary edge was encountered. Insert (ctrx,ctry) into the * triangulation (whether or not kt is contained in R). It may have * to be deleted if it encroaches upon a boundary edge. */ if (nv+1 > nvmax) trealloc(nv+1); insertv(kt, b1, b2, b3); nsp++; /* * Apply a sequence of edge swaps to restore the Delaunay circle * property. The edge indices are stored in lste. */ neq0 = neq; update(nv-1, 0, &ns, lste); if (neq == neq0) { /* * No boundary edges were encroached. Repeat the loop on triangles * containing nv-1, this time with bad triangles of R added to Qt. */ update(nv-1, 1, &ns, lste); } else { /* * One or more boundary edges were encroached upon by (ctrx,ctry). * Their indices have been enqueued. The swaps are reversed, * and (ctrx,ctry) is deleted. */ for (k = ns-1; k >= 0; k--) { unswap(lste[k]); } if (!delvtx(nv-1)) { printf("glmesh2: Failure to reverse a vertex " "insertion in refineRd.\n"); exit(1); } nsp--; } } else { /* * Enqueue edge ke. */ kt = ledg[ke]; if (ltri[kt][6] == ke) i = 0; else if (ltri[kt][7] == ke) i = 1; else i = 2; kv1 = ltri[kt][(i+1)%3]; kv2 = ltri[kt][(i+2)%3]; if (neq == neqmax) reallocQe(); qeback = (qeback+1)%neqmax; Qe[qeback].index_e = ke; Qe[qeback].index_v1 = kv1; Qe[qeback].index_v2 = kv2; neq++; } if (ke >= 0 || neq > neq0) { /* * The off-center (ctrx,ctry) was not inserted because it encroached * upon a boundary edge. Put the triangle back in the priority queue * for another try, and split the encroached edges. */ insertQt(kt0, imin, emin); splitEdges(1); /* * Test for failure to converge. */ if (neq) { printf("\n splitEdges aborted with neq = %d\n",neq); break; } } } /* * Test for failure to converge. */ if (neq || ntq) { printf("\n refineRd aborted with neq = %d, ntq = %d\n", neq, ntq); } /* * Free storage. */ free(Qt); free(Qtp); free(Qe); Qt = NULL; Qtp = NULL; Qe = NULL; return; } /* End of refineRd */ /* *-------------------------------------------------------------------- * * refineRe: Refine R by adding a new vertex at the midpoint of each * edge in R. The new vertex is connected to the end- * points of the original edge and to the one or two * vertices opposite the original edge, thus replacing * each triangle involved by two new triangles. Each * triangle contained in R is partitioned into four equal- * area subtriangles by three interior edges. Following * the additions, an edge swap is applied to the first * interior edge that was added to each triangle in R so * that, if the original triangle was equilateral, the * four new triangles, which are all similar to the * original triangle, are also equilateral. * * glmesh2 functions called: initR, split, swap, trealloc * *-------------------------------------------------------------------- */ static void refineRe(void) { int i, k, ke, ken, kt, ktn, kv0, kv3, ne0, ner, ner0, nv0; /* * Store listR if necessary. */ if (newR) initR(); /* * Store a list, lste, of the indices of the edges contained in R. * The count is stored in ner. */ ner = 0; for (ke = 0; ke < ne; ke++) { kt = ledg[ke]; if (ltri[kt][6] == ke) i = 1; else if (ltri[kt][7] == ke) i = 2; else i = 0; if (!listR[ltri[kt][i]] || !listR[ltri[kt][(i+1)%3]]) continue; lste[ner] = ke; ner++; } if (!ner) return; /* * Push nv onto rstack for Function unrefine. */ if (nref == nrefmax) { nrefmax *= 2; rstack = realloc(rstack, nrefmax*sizeof *rstack); if (rstack == NULL) { printf("glmesh2: Unable to allocate sufficient memory " "for rstack.\n"); exit(1); } } rstack[nref] = nv; nref++; if (nv + ner > nvmax) { /* * Increase nvmax, and re-allocate storage for vtxy, uv, lvtx, ltri, * ledg, lnrv, lcrv, lste, listn, and listR. */ trealloc(nv+ner); } /* * Initialization: * * ne0 = Initial edge count, * nv0 = Initial vertex count. * ner0 = Initial value of ner (number of edges in R). */ ne0 = ne; nv0 = nv; ner0 = ner; /* * Loop on edges ke in R. */ for (k = 0; k < ner0; k++) { ke = lste[k]; split(ke, 0.5, 0.5); /* * The first edges added to the interior of the triangles that * were split have endpoint indices nv-1 and kv0 < nv0. If either * triangle is in R, the corresponding edge index is appended to * lste, indicating that it is to be swapped. */ kt = ledg[ke]; if (ltri[kt][6] == ke) i = 0; else if (ltri[kt][7] == ke) i = 1; else i = 2; kv0 = ltri[kt][i]; ktn = ltri[kt][i+3]; if (listR[kv0] && kv0 < nv0) { ken = ne-1; if (ktn >= 0) ken--; lste[ner] = ken; ner++; } /* * The first edge added to the interior of triangle ktn, if ktn is * in R, must also be added to the list. */ if (ktn >= 0) { if (ltri[ktn][3] == kt) i = 0; else if (ltri[ktn][4] == kt) i = 1; else i = 2; kv3 = ltri[ktn][i]; if (listR[kv3] && kv3 < nv0) { lste[ner] = ne-1; ner++; } } } /* end of loop on edges */ /* * Apply swaps to the edges that were added to lste. */ for (k = ner0; k < ner; k++) { ke = lste[k]; swap(ke); } return; } /* End of refineRe */ /* *-------------------------------------------------------------------- * * refineRt: Refine R by adding triangle barycenters. Refer to * Function insertv. * * glmesh2 functions called: initR, insertv, trealloc * *-------------------------------------------------------------------- */ static void refineRt(void) { double b = 1.0/3.0; /* Barycentric coordinates of new vertex */ int k, kt, ntr; /* * Store listR if necessary. */ if (newR) initR(); /* * Store a list, lste, of triangles contained in R. * The count is stored in ntr. */ ntr = 0; for (kt = 0; kt < nt; kt++) { if (!listR[ltri[kt][0]] || !listR[ltri[kt][1]] || !listR[ltri[kt][2]]) continue; lste[ntr] = kt; ntr++; } if (!ntr) return; /* * Push nv onto rstack for Function unrefine. */ if (nref == nrefmax) { nrefmax *= 2; rstack = realloc(rstack, nrefmax*sizeof *rstack); if (rstack == NULL) { printf("glmesh2: Unable to allocate sufficient memory " "for rstack.\n"); exit(1); } } rstack[nref] = nv; nref++; if (nv + ntr > nvmax) { /* * Increase nvmax, and re-allocate storage for vtxy, uv, lvtx, ltri, * ledg, lnrv, lcrv, lste, listn, and listR. */ trealloc(nv+ntr); } /* * Loop on triangles in R, inserting the barycenters. */ for (k = 0; k < ntr; k++) { insertv(lste[k], b, b, b); } return; } /* End of refineRt */ /* *-------------------------------------------------------------------- * * removeRt: Remove the triangles in R from the triangle mesh (and * from R), creating new boundary curves if necessary. It * may not be possible to remove all triangles if removal * would render the mesh invalid. * * glmesh2 functions called: deltri, initR * *-------------------------------------------------------------------- */ static void removeRt(void) { int btest = 1; /* Flag specifying a test for a boundary edge */ int k; /* Index for lste */ int kt; /* Triangle index */ int nd = 1; /* Flag specifying another sweep required */ int ntr = 0; /* Number of triangles in R initially */ int nv0 = nv; /* Initial vertex count */ /* * Store listR if necessary. */ if (newR) initR(); /* * Store a list, lste, of triangles contained in R. * The count is stored in ntr. */ for (kt = 0; kt < nt; kt++) { if (!listR[ltri[kt][0]] || !listR[ltri[kt][1]] || !listR[ltri[kt][2]]) continue; lste[ntr] = kt; ntr++; } if (!ntr) return; /* * Outer loop on sweeps. */ while (nd) { nd = 0; k = 0; /* * Inner loop on triangles in R. */ while (k < ntr) { kt = lste[k]; if (!btest || (ltri[kt][3] < 0 || ltri[kt][4] < 0 || ltri[kt][5] < 0)) { if (deltri(kt)) { nd++; if (!btest) btest = 1; } else { k++; } } else { k++; } } if (!nd && btest) { btest = 0; nd = 1; } } /* * The refinement stack is emptied if a vertex was deleted. */ if (nv < nv0) nref = 0; return; } /* End of removeRt */ /* *-------------------------------------------------------------------- * * reshape: Set up a viewport, construct a projection, and * initialize the modelview matrix. * * Global variables required: aspect, dx, dy, sx, sy, xc, yc * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void reshape(int w, int h) { GLdouble r; /* Aspect ratio of the view volume */ int imax, imin, jmax, jmin; /* Window coordinates of viewport */ int ni, nj; /* Width and height of viewport */ if (aspect) { /* * Compute viewport corner coordinates such that the aspect ratio r * of the viewing volume is preserved (assuming square pixels) and * the viewport is centered in the window and as large as possible. */ r = dx/dy; if (r <= (GLdouble) w/(GLdouble) h) { imin = (int) ((w-r*h)/2.0); imax = imin + (int) (r*h+0.5); if (imax > w) imax = w; jmin = 0; jmax = h; } else { imin = 0; imax = w; jmin = (int) ((h-w/r)/2.0); jmax = jmin + (int) (w/r+0.5); } } else { imin = 0; imax = w; jmin = 0; jmax = h; } ni = imax-imin; nj = jmax-jmin; /* * Store global variable sx and sy, and create the viewport mapping. */ sx = dx/(GLdouble) ni; sy = dy/(GLdouble) nj; glViewport(imin, jmin, ni, nj); glMatrixMode(GL_PROJECTION); glLoadIdentity(); /* * Orthographic projection. */ gluOrtho2D(xc-dx/2.0, xc+dx/2.0, yc-dy/2.0, yc+dy/2.0); glMatrixMode(GL_MODELVIEW); glLoadIdentity(); return; } /* End of reshape */ /* *-------------------------------------------------------------------- * * specialKey: Move center (xc,yc) in response to arrow keys, zoom * in or out on Page Up or Page Down keys, or restore * default viewing volume if Home key is pressed. * * Global variables required: dx, dy, xc, yc, xmin, ymin, xmax, ymax * * glmesh2 function called: reshape * *-------------------------------------------------------------------- */ static void specialKey(int key, int x, int y) { switch (key) { case GLUT_KEY_UP: if (glutGetModifiers() == GLUT_ACTIVE_SHIFT) { yc += 0.01*(ymax-ymin); } else { yc += 0.05*(ymax-ymin); } break; case GLUT_KEY_DOWN: if (glutGetModifiers() == GLUT_ACTIVE_SHIFT) { yc -= 0.01*(ymax-ymin); } else { yc -= 0.05*(ymax-ymin); } break; case GLUT_KEY_RIGHT: if (glutGetModifiers() == GLUT_ACTIVE_SHIFT) { xc += 0.01*(xmax-xmin); } else { xc += 0.05*(xmax-xmin); } break; case GLUT_KEY_LEFT: if (glutGetModifiers() == GLUT_ACTIVE_SHIFT) { xc -= 0.01*(xmax-xmin); } else { xc -= 0.05*(xmax-xmin); } break; case GLUT_KEY_PAGE_UP: /* Zoom in */ dx *= 0.8; dy *= 0.8; break; case GLUT_KEY_PAGE_DOWN: /* Zoom out */ dx *= 1.25; dy *= 1.25; break; case GLUT_KEY_HOME: /* Restore default viewing volume */ dx = 1.2*(xmax-xmin); dy = 1.2*(ymax-ymin); xc = (xmin+xmax)/2.0; yc = (ymin+ymax)/2.0; break; default: break; } reshape(glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT)); glutPostRedisplay(); return; } /* End of specialKey */ /* *-------------------------------------------------------------------- * * split: Add a new vertex in the interior of an edge, splitting the * edge into a pair of edges, and replacing each triangle * containing the edge by two triangles. The uv values at * the new vertex are obtained by linear interpolation. * * On input: ke = index of the edge to be split. * * b1 and b2 are the local coordinates of the new vertex * with respect to the edge. It is assumed without a * test that b1 and b2 are nonnegative and add to 1. * * It is necessary that nv < nvmax so that there is sufficient room * in arrays vtxy, ltri, lvtx, and ledg to add another vertex, two * more triangles, and three more edges. No variables are altered * if nv >= nvmax. * * The new vertex is flagged as being contained in R (listR[nv-1] = * 1) iff both endpoint vertices of edge ke are in R. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void split(int ke, double b1, double b2) { int i, j, ke1, ke2 = 0, kt, kt1, kt2, ktn, kv0, kv1, kv2, kv3 = 0; if (nv >= nvmax) return; /* * Set (kv0,kv1,kv2) to the vertex indices of the triangle kt that * contains ke, where kv0 = ltri[kt][i] and ke = ltri[kt][i+6]. */ kt = ledg[ke]; if (ltri[kt][6] == ke) i = 0; else if (ltri[kt][7] == ke) i = 1; else i = 2; kv0 = ltri[kt][i]; kv1 = ltri[kt][(i+1)%3]; kv2 = ltri[kt][(i+2)%3]; /* * Insert the new vertex by appending its Cartesian coordinates * and uv values. */ vtxy[nv][0] = b1*vtxy[kv1][0] + b2*vtxy[kv2][0]; vtxy[nv][1] = b1*vtxy[kv1][1] + b2*vtxy[kv2][1]; uv[nv][0] = b1*uv[kv1][0] + b2*uv[kv2][0]; uv[nv][1] = b1*uv[kv1][1] + b2*uv[kv2][1]; /* * Set ktn to the index of the triangle opposite kv0, and set kt1 * and ke1 to the indices of the triangle and edge opposite kv1. */ ktn = ltri[kt][i+3]; kt1 = ltri[kt][(i+1)%3+3]; ke1 = ltri[kt][(i+1)%3+6]; /* * If i = j = 0, the initial ltri rows are * * kt: (kv0, kv1, kv2, ktn, kt1, x, ke, ke1, x) * ktn: (kv3, kv2, kv1, kt, x, kt2, ke, x, ke2) if ktn > 0, * * and the new columns, following the split are * * kt: (kv0, kv1, nv, ktn, nt, x, ke, ne+1, x) * nt: (kv0, nv, kv2, nt+1, kt1, kt, ne, ke1, ne+1) * nt+1: (kv3, kv2, nv, nt, ktn, kt2, ne, ne+2, ke2) * ktn: (kv3, nv, kv1, kt, x, nt+1, ke, x, ne+2) * * where the entries labeled x retain their values and need not be * referenced, nt+1 = -1 if ktn = -1, and the last two rows only * exist if ktn >= 0. The changes implicit in the above lists are * identical for other cyclical orderings (other values of i and j). */ ltri[kt][(i+2)%3] = nv; ltri[kt][(i+1)%3+3] = nt; ltri[kt][(i+1)%3+6] = ne+1; ltri[nt][0] = kv0; ltri[nt][1] = nv; ltri[nt][2] = kv2; ltri[nt][3] = -1; ltri[nt][4] = kt1; ltri[nt][5] = kt; ltri[nt][6] = ne; ltri[nt][7] = ke1; ltri[nt][8] = ne+1; kt2 = -1; if (ktn >= 0) { /* * Find j such such that kv3 = ltri[ktn][j] is the vertex opposite * edge ke, and set kt2 and ke2 to the triangle and edge opposite * kv1 in ktn. */ if (ltri[ktn][0] == kv1) j = 1; else if (ltri[ktn][1] == kv1) j = 2; else j = 0; kv3 = ltri[ktn][j]; kt2 = ltri[ktn][(j+2)%3+3]; ke2 = ltri[ktn][(j+2)%3+6]; ltri[nt][3] = nt+1; ltri[ktn][(j+1)%3] = nv; ltri[ktn][(j+2)%3+3] = nt+1; ltri[ktn][(j+2)%3+6] = ne+2; ltri[nt+1][0] = kv3; ltri[nt+1][1] = kv2; ltri[nt+1][2] = nv; ltri[nt+1][3] = nt; ltri[nt+1][4] = ktn; ltri[nt+1][5] = kt2; ltri[nt+1][6] = ne; ltri[nt+1][7] = ne+2; ltri[nt+1][8] = ke2; } /* * Correct neighboring triangle ltri entries. */ if (kt1 >= 0) { if (ltri[kt1][3] == kt) ltri[kt1][3] = nt; else if (ltri[kt1][4] == kt) ltri[kt1][4] = nt; else ltri[kt1][5] = nt; } if (kt2 >= 0) { if (ltri[kt2][3] == ktn) ltri[kt2][3] = nt+1; else if (ltri[kt2][4] == ktn) ltri[kt2][4] = nt+1; else ltri[kt2][5] = nt+1; } /* * Update lvtx. */ if (lvtx[kv2] == kt) lvtx[kv2] = nt; if (lvtx[kv2] == ktn) lvtx[kv2] = nt+1; if (ktn >= 0) { if (lvtx[kv3] == ktn) lvtx[kv3] = nt+1; } lvtx[nv] = nt; /* * Update ledg, preserving the property that the larger of the two * triangle indices is used. */ ledg[ne] = nt; ledg[ne+1] = nt; ledg[ke1] = nt; if (ktn >= 0) { ledg[ne] = nt+1; ledg[ne+2] = nt+1; ledg[ke2] = nt+1; } /* * Update the counts. */ nv++; if (ktn < 0) { nt++; ne += 2; nb++; } else { nt += 2; ne += 3; } /* * Update listn and listR. */ listn[nv-1] = 0; listR[nv-1] = listR[kv1] && listR[kv2]; return; } /* End of split */ /* *-------------------------------------------------------------------- * * splitEdges: Split all encroached boundary edges in R by inserting * their midpoints (or points close to the midpoints) * into the triangle mesh: used by Function refineRd. * * On input: The boundary edges that are encroached upon by their * opposite vertices must be stored in Qe. * * tritest is a flag with value 1 iff newly created * triangles are to be tested for inclusion in Qt by * function badTriangle. * * glmesh2 functions called: encroached, split, trealloc, update * *-------------------------------------------------------------------- */ static void splitEdges(unsigned char tritest) { double b1, b2, d, d0, dx, dy; int i, ke, kt, kv1, kv2, ns; /* * Loop on edges in Qe. */ while (neq && nsp < nspmax) { /* * Dequeue edge ke with endpoint indices kv1 and kv2. */ ke = Qe[qefront].index_e; kv1 = Qe[qefront].index_v1; kv2 = Qe[qefront].index_v2; qefront = (qefront+1)%neqmax; neq--; /* * Bypass edges that have been removed from the triangle mesh. */ kt = ledg[ke]; if (ltri[kt][6] == ke) i = 0; else if (ltri[kt][7] == ke) i = 1; else i = 2; if (ltri[kt][(i+1)%3] != kv1 || ltri[kt][(i+2)%3] != kv2) continue; /* * Reallocate storage for the triangle mesh if necessary. */ if (nv+1 > nvmax) trealloc(nv+1); /* * Split edge ke. If the triangle containing ke has two boundary * edges, then the midpoint of ke could encroach upon the other * boundary edge, causing it to be split with a new vertex that * again encroached upon a boundary edge (one of the two halves * of ke), ad infinitum. Rather than splitting at the midpoint in * this case, the splitting point is taken to be the intersection * of the edge with one of the concentric circles, centered at the * vertex shared by the two boundary edges, and having radius equal * to a power of two. By choosing the nearest power of 2 to half * the edge length, the ratio of lengths of the two edge portions * is at most 2 to 1: the local coordinates of the split point are * at least 1/3. This scheme precludes the possibility of an * infinite number of splits. */ if (ltri[kt][(i+1)%3+3] < 0 || ltri[kt][(i+2)%3+3] < 0) { dx = vtxy[kv2][0] - vtxy[kv1][0]; dy = vtxy[kv2][1] - vtxy[kv1][1]; d0 = 0.5*sqrt(dx*dx + dy*dy); d = pow(2.0, floor(log(d0)/log(2.0))); if (d0 > 1.5*d) d = 2.0*d; b1 = 0.5*d/d0; b2 = 1.0 - b1; if (ltri[kt][(i+2)%3+3] < 0) { b2 = b1; b1 = 1.0 - b2; } } else { b1 = 0.5; b2 = 0.5; } /* * Insert the splitting point of edge ke into the triangle mesh as * vertex nv-1, and and update R. The two portions of edge ke are * assigned edge indices ke and ne-2. */ split(ke, b1, b2); nsp++; /* * Apply the necessary edge swaps to restore the empty circumcircle * property, and update Qe for the edges opposite the new vertex, and * update Qt if tritest = 1. The swapped edge indices are stored in * lste but are not used. */ update(nv-1, tritest, &ns, lste); /* * Test the two new boundary edges incident on vertex nv-1 for * encroachment. */ if (encroached(ke, &kv1, &kv2)) { if (neq == neqmax) reallocQe(); qeback = (qeback+1)%neqmax; Qe[qeback].index_e = ke; Qe[qeback].index_v1 = kv1; Qe[qeback].index_v2 = kv2; neq++; } ke = ne-2; if (encroached(ke, &kv1, &kv2)) { if (neq == neqmax) reallocQe(); qeback = (qeback+1)%neqmax; Qe[qeback].index_e = ke; Qe[qeback].index_v1 = kv1; Qe[qeback].index_v2 = kv2; neq++; } } return; } /* End of splitEdges */ /* *-------------------------------------------------------------------- * * storeV: Compute and store global parameter values that depend on * the input data: dx, dy, xc, xmax, xmin, yc, ymax, and * ymin, and allocate storage for R and rstack. The initial * values of scamax and offscale are also computed. * * Global variables required: * * min_angle = Minimum angle defining a bad triangle for Delaunay * refinement (Function refineRd). * * npmax = Initial size for R. * * nrefmax = Initial size for rstack. * * nv = Number of vertices in the triangle mesh. * * R = Array dimensioned npmax by 2 containing the Cartesian * coordinates of a sequence of vertices defining the * boundary of the polygonal region. * * rstack = Array of length nrefmax containing a sequence of * vertex counts nv associated with mesh refinements. * * vtxy = Array dimensioned nv by 2 containing the Cartesian * coordinates of the triangle mesh vertices. * * glmesh2 function called: bbox * *------------------------------------------------------------------- */ static void storeV(void) { double c; /* * Allocate storage. */ R = malloc(npmax*sizeof *R); rstack = malloc(nrefmax*sizeof *rstack); /* * Compute the range of vertex coordinates (bounding box corner * coordinates): [xmin,xmax] X [ymin,ymax]. */ bbox(); /* * Store the viewing volume dimensions and center. */ dx = 1.2*(xmax-xmin); dy = 1.2*(ymax-ymin); xc = (xmin+xmax)/2.0; yc = (ymin+ymax)/2.0; /* * Compute initial values that depend on min_angle. */ c = cos(acos(-1.0)*min_angle/180.0); scamax = c*c; offscale = 0.48*sqrt((1.0+c)/(1.0-c)); return; } /* End of storeV */ /* *-------------------------------------------------------------------- * * swap: Modify a triangle mesh by swapping diagonal edges in a con- * vex quadrilateral defined by a pair of adjacent triangles * (triangles that share a common edge). For a given set of * vertices, any triangle mesh can be transformed into any * other triangle mesh by a sequence of such swaps. * * On input: * * ke = Index of the edge to be swapped. * * It is assumed that the quadrilateral is convex. No variables are * altered, and the return value is 0 if edge ke is not shared by two * triangles, or if the vertices opposite edge ke are already * adjacent. Otherwise the return value is 1. * * Global variables required: ledg, ltri, lvtx * * glmesh2 function called: adjacent * *-------------------------------------------------------------------- */ static int swap(int ke) { int i, j, ken1, ken2, kt1, kt2, ktn1, ktn2, kv1, kv2, kvn1, kvn2; /* * Set kt1 and kt2 to the adjacent triangles with shared edge ke, * set (kv1,kvn1,kvn2) and (kv2,kvn2,kvn1) to the vertex indices of * triangles kt1 and kt2, respectively, and set i and j to the column * indices of kv1 in kt1 and kv2 in kt2, respectively. */ kt1 = ledg[ke]; if (ltri[kt1][6] == ke) i = 0; else if (ltri[kt1][7] == ke) i = 1; else i = 2; kt2 = ltri[kt1][i+3]; /* * Test for a boundary edge. */ if (kt2 < 0) return 0; if (ltri[kt2][6] == ke) j = 0; else if (ltri[kt2][7] == ke) j = 1; else j = 2; kv1 = ltri[kt1][i]; kv2 = ltri[kt2][j]; kvn1 = ltri[kt1][(i+1)%3]; kvn2 = ltri[kt2][(j+1)%3]; ktn1 = ltri[kt2][(j+1)%3+3]; ktn2 = ltri[kt1][(i+1)%3+3]; ken1 = ltri[kt2][(j+1)%3+6]; ken2 = ltri[kt1][(i+1)%3+6]; /* * If i = j = 0, the initial ltri rows are * * kt1 = (kv1,kvn1,kvn2,kt2,ktn2,x,ke,ken2,x) * kt2 = (kv2,kvn2,kvn1,kt1,ktn1,x,ke,ken1,x) * * and the new rows, following the swap are * * kt1 = (kv1,kvn1,kv2,ktn1,kt2,x,ken1,ke,x) * kt2 = (kv2,kvn2,kv1,ktn2,kt1,x,ken2,ke,x) * * where the entries labeled x retain their values and need not be * referenced. * * Test for kv1 already adjacent to kv2. */ if (adjacent(kv1,kv2)) { printf("swap: WARNING: Failed attempt to connect vertices " "that are already adjacent.\n"); return 0; } /* * Change ltri entries. */ ltri[kt1][(i+2)%3] = kv2; ltri[kt2][(j+2)%3] = kv1; ltri[kt1][i+3] = ktn1; ltri[kt2][j+3] = ktn2; ltri[kt1][(i+1)%3+3] = kt2; ltri[kt2][(j+1)%3+3] = kt1; ltri[kt1][i+6] = ken1; ltri[kt2][j+6] = ken2; ltri[kt1][(i+1)%3+6] = ke; ltri[kt2][(j+1)%3+6] = ke; /* * Correct neighboring triangle ltri entries if necessary. */ if (ktn1 >= 0) { if (ltri[ktn1][3] == kt2) ltri[ktn1][3] = kt1; else if (ltri[ktn1][4] == kt2) ltri[ktn1][4] = kt1; else ltri[ktn1][5] = kt1; } if (ktn2 >= 0) { if (ltri[ktn2][3] == kt1) ltri[ktn2][3] = kt2; else if (ltri[ktn2][4] == kt1) ltri[ktn2][4] = kt2; else ltri[ktn2][5] = kt2; } /* * Adjust lvtx entries. */ if (lvtx[kvn1] == kt2) lvtx[kvn1] = kt1; if (lvtx[kvn2] == kt1) lvtx[kvn2] = kt2; /* * Adjust ledg entries, preserving the property that the larger of * the two triangle indices is used. */ ledg[ken1] = kt1; if (ktn1 > kt1) ledg[ken1] = ktn1; ledg[ken2] = kt2; if (ktn2 > kt2) ledg[ken2] = ktn2; return 1; } /* End of swap */ /* *-------------------------------------------------------------------- * * swapEdges: Mouse callback (callback triggered by a mouse button * press or release with the mouse position in the * primary window) used to select an edge to be swapped. * * This function swaps the edge selected by a left button press if * the mouse position is close enough to an edge, and the edge can * be swapped. * * glmesh2 functions called: plDistance, swap, xlate * *-------------------------------------------------------------------- */ static void swapEdges(int button, int state, int x, int y) { GLdouble wx, wy, wz; /* Computed object coordinates */ GLboolean intr; GLdouble tols; GLint i, ke, kt, ktn, kv0, kv1 = 0, kv2 = 0, kv3; int swp; if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) { xlate(x,y,0.0, &wx,&wy,&wz); /* * Compute a tolerance tols on the squared distance from the * mouse position (wx,wy) to the edge. */ tols = (sx <= sy) ? sx : sy; /* tols = min(sx, sy) */ tols = 16.0*tols*tols; /* 4-pixel tolerance */ /* * Loop on edges kv1-kv2. */ intr = 0; for (kt = 0; kt < nt; kt++) { for (i = 0; i < 3; i++) { ktn = ltri[kt][i+3]; if (ktn > kt) continue; kv1 = ltri[kt][(i+1)%3]; kv2 = ltri[kt][(i+2)%3]; if ((intr = plDistance(wx, wy, vtxy[kv1][0], vtxy[kv1][1], vtxy[kv2][0], vtxy[kv2][1], tols))) break; } if (intr) break; } if (!intr) return; if (ktn < 0) { printf("\n A boundary edge cannot be swapped.\n"); return; } kv0 = ltri[kt][i]; ke = ltri[kt][i+6]; if (ltri[ktn][6] == ke) kv3 = ltri[ktn][0]; else if (ltri[ktn][7] == ke) kv3 = ltri[ktn][1]; else kv3 = ltri[ktn][2]; /* * Edge ke can be swapped iff kv1 strictly Left kv3 -> kv0 and * kv2 strictly Left kv0 -> kv3. */ swp = ((vtxy[kv0][0]-vtxy[kv3][0])*(vtxy[kv1][1]-vtxy[kv3][1]) > (vtxy[kv0][1]-vtxy[kv3][1])*(vtxy[kv1][0]-vtxy[kv3][0]) && (vtxy[kv3][0]-vtxy[kv0][0])*(vtxy[kv2][1]-vtxy[kv0][1]) > (vtxy[kv3][1]-vtxy[kv0][1])*(vtxy[kv2][0]-vtxy[kv0][0])); if (swp) { if (swap(ke)) { glutPostRedisplay(); } else { swp = 0; } } if (!swp) printf("\n That edge cannot be swapped.\n"); } return; } /* End of swapEdges */ /* *-------------------------------------------------------------------- * * swptst: Delaunay circumcircle test. * * This function applies the circumcircle test to a quadrilateral * defined by a pair of adjacent triangles. The diagonal edge * (shared triangle side) should be swapped for the other diagonl if * and only if the fourth vertex is strictly interior to the circum- * circle of one of the triangles (the decision is independent of the * choice of triangle). Equivalently, the diagonal is chosen to * maximize the smallest of the six interior angles over the two * pairs of possible triangles (and the decision is for no swap if * the quadrilateral is not strictly convex). * * When the four vertices are nearly cocircular (the neutral case), * the preferred decision is no swap -- in order to avoid unneces- * sary swaps and, more important, to avoid cycling in a sequence of * sweeps (through a set of triangles) terminated by the global * empty circumcircle property (no circumcircle of a triangle in * the set contains a vertex in its interior). Thus, a tolerance * swtol is used to define 'nearness' to the neutral case. * * On input: in1, in2, io1, and io2 are the vertex indices (indices * for vtxy) of the quadrilateral with triangles defined * by CCW-ordered triples (io1,io2,in1) and (io2,io1,in2) * so that io1-io2 is the shared side to be replaced by * edge in1-in2. * * On output: the return value is 1 iff the shared edge io1-io2 * should be swapped for the other diagonal in1-in2. * * Global variable required: vtxy * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static unsigned char swptst(int in1, int in2, int io1, int io2) { double cos1, cos2, dx11, dx12, dx21, dx22, dy11, dy12, dy21, dy22, sin1, sin2, sin12; double swtol = 2.2e-15; /* 10*(Machine precision) */ /* * Compute components of the directed edges anchored at in1 and in2. */ dx11 = vtxy[io1][0] - vtxy[in1][0]; dx12 = vtxy[io2][0] - vtxy[in1][0]; dx21 = vtxy[io1][0] - vtxy[in2][0]; dx22 = vtxy[io2][0] - vtxy[in2][0]; dy11 = vtxy[io1][1] - vtxy[in1][1]; dy12 = vtxy[io2][1] - vtxy[in1][1]; dy21 = vtxy[io1][1] - vtxy[in2][1]; dy22 = vtxy[io2][1] - vtxy[in2][1]; /* * Compute inner products, proportional cos(A1) and cos(A2), where * A1 and A2 are the angles at in1 and in2, respectively. */ cos1 = dx11*dx12 + dy11*dy12; cos2 = dx21*dx22 + dy21*dy22; /* * The diagonals should be swapped iff (A1+A2) > 180 degrees. The * following two tests ensure numerical stability: the decision must * be FALSE when both angles are close to 0, and TRUE when both * angles are close to 180 degrees. */ if (cos1 >= 0.0 && cos2 >= 0.0) return 0; if (cos1 < 0.0 && cos2 < 0.0) return 1; /* * Compute vector cross products (Z-components) proportional to * sin(A1) and sin(A2), and sin12 proportional to sin(A1+A2). */ sin1 = dx11*dy12 - dx12*dy11; sin2 = dx22*dy21 - dx21*dy22; sin12 = sin1*cos2 + cos1*sin2; return (sin12 < -swtol); } /* End of swptst */ /* *-------------------------------------------------------------------- * * trealloc: Re-allocate storage for the triangle mesh. * * On input: nvnew is the required (minimum) value of nvmax. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void trealloc(int nvnew) { /* * Increase nvmax, and re-allocate storage for vtxy, uv, lvtx, ltri, * ledg, lnrv, lcrv, lste, listn, and listR. */ nvmax = 1.1*(nvnew); vtxy = realloc(vtxy, nvmax*sizeof *vtxy); uv = realloc(uv, nvmax*sizeof *uv); lvtx = realloc(lvtx, nvmax*sizeof *lvtx); ltri = realloc(ltri, 2*nvmax*sizeof *ltri); ledg = realloc(ledg, 3*nvmax*sizeof *ledg); lnrv = realloc(lnrv, nvmax*sizeof *lnrv); lcrv = realloc(lcrv, (nvmax/3)*sizeof *lcrv); lste = realloc(lste, 3*nvmax*sizeof *lste); listn = realloc(listn, nvmax*sizeof *listn); listR = realloc(listR, nvmax*sizeof *listR); if (vtxy == NULL || uv == NULL || lvtx == NULL || ltri == NULL || ledg == NULL || lnrv == NULL || lcrv == NULL || lste == NULL || listn == NULL || listR == NULL) { printf("glmesh2: Unable to allocate sufficient memory.\n"); exit(1); } return; } /* End of trealloc */ /* *-------------------------------------------------------------------- * * trfind: Locate a point relative to the triangle mesh. * * Note that, rather than solving the general point-location problem, * which is difficult in a non-convex mesh, this function is designed * for the case that the target point p is in the circumcircle of the * starting triangle kt0. In particular, it is assumed that at most * one of the barycentric coordinates of p with respect to triangle * kt0 is negative. * * On input: kt0 indexes a triangle in which the search begins, and * (xp,yp) is the target point. * * On output: If a line segment joining p to an interior point of * triangle kt0 intersects a boundary edge, kep indexes * that boundary edge (or the one closest to kt0 if there * is more than one), the return value is 0, and no other * output parameters are altered. Otherwise, the return * value is positive, ktp indexes a triangle containing * p, the barycentric coordinates of p with respect to * ktp are stored in b0p, b1p, and b2p (in the same order * as the vertices), and kep is unaltered. If triangle * ktp is null, p is taken to be the barycenter, a * warning message is printed, and the return value is 2. * Otherwise, the return value is 1. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static int trfind(int kt0, double xp, double yp, int *kep, int *ktp, double *b0p, double *b1p, double *b2p) { double b, b0, b1, b2, x0, y0; int i, kt, ktn, kv0, kv1, kv2, kv3, kv4; /* * Initialize kt = kt0 with vertices (kv0,kv1,kv2). */ kt = kt0; kv0 = ltri[kt][0]; kv1 = ltri[kt][1]; kv2 = ltri[kt][2]; /* * Test for p in kt0 (found = 1), or relabel the vertex indices * so that v0-p intersects v1-v2, where v0, v1, and v2 are the * vertices indexed by kv0, kv1, and kv2. */ i = 0; b0 = (vtxy[kv1][0]-xp)*(vtxy[kv2][1]-yp) - (vtxy[kv2][0]-xp)*(vtxy[kv1][1]-yp); if (b0 >= 0.0) { b1 = (vtxy[kv2][0]-xp)*(vtxy[kv0][1]-yp) - (vtxy[kv0][0]-xp)*(vtxy[kv2][1]-yp); if (b1 < 0) { i = 1; kv1 = kv2; kv2 = kv0; kv0 = ltri[kt][i]; b0 = b1; } else { b2 = (vtxy[kv0][0]-xp)*(vtxy[kv1][1]-yp) - (vtxy[kv1][0]-xp)*(vtxy[kv0][1]-yp); if (b2 < 0) { i = 2; kv2 = kv1; kv1 = kv0; kv0 = ltri[kt][i]; b0 = b2; } } } /* * Initialization for edge-hopping loop with the following loop * invariant: v0-p intersects v1-v2, where triangle kt has * vertices v1, v2, and v3 indexed by kv1, kv2, and kv3 = ltri[kt][i]. */ x0 = vtxy[kv0][0]; y0 = vtxy[kv0][1]; kv3 = kv0; while (b0 < 0.0) { ktn = ltri[kt][i+3]; if (ktn < 0) { /* * v1-v2 is a boundary edge. */ *kep = ltri[kt][i+6]; return 0; } /* * Advance to the next edge. Neighboring triangle ktn has vertex * indices (kv2,kv1,kv4). */ if (ltri[ktn][3] == kt) i = 0; else if (ltri[ktn][4] == kt) i = 1; else i = 2; kt = ktn; kv4 = ltri[kt][i]; /* * Test for p Left v0 -> v4. */ if ((vtxy[kv4][0]-x0)*(yp-y0) >= (xp-x0)*(vtxy[kv4][1]-y0)) { kv3 = kv1; kv1 = kv4; i = (i+2)%3; } else { kv3 = kv2; kv2 = kv4; i = (i+1)%3; } /* * Compute b0 = (v2-v1) X (p-v1) = (v1-p) X (v2-p) >= 0 iff * p Left v1 -> v2. */ b0 = (vtxy[kv1][0]-xp)*(vtxy[kv2][1]-yp) - (vtxy[kv2][0]-xp)*(vtxy[kv1][1]-yp); } /* * p is in triangle kt. Compute the barycentric coordinates of p * with respect to the triangle. */ b1 = (vtxy[kv2][0]-xp)*(vtxy[kv3][1]-yp) - (vtxy[kv3][0]-xp)*(vtxy[kv2][1]-yp); b2 = (vtxy[kv3][0]-xp)*(vtxy[kv1][1]-yp) - (vtxy[kv1][0]-xp)*(vtxy[kv3][1]-yp); if (i == 1) { b = b0; b0 = b2; b2 = b1; b1 = b; } else if (i == 2) { b = b0; b0 = b1; b1 = b2; b2 = b; } b = b0+b1+b2; if (b != 0.0) { b = 1.0/b; *ktp = kt; *b0p = b0*b; *b1p = b1*b; *b2p = b2*b; return 1; } else { /* * b0+b1+b2 = 0. Return the barycenter of triangle kt. */ printf("\n trfind: insertion at barycenter of null triangle " "%d.\n", kt); *ktp = kt; b0 = 1.0/3.0; *b0p = b0; *b1p = b0; *b2p = b0; return 2; } } /* End of trfind */ /* *-------------------------------------------------------------------- * * trmtst: Triangulation data structure integrity test. * * This function tests for the following properties: * * 1) Counts nc, nb, nv, nn, nt, and ne are consistent. * 2) Vertex indices in ltri, lcrv, and lnrv are in the range 0 * to nv-1. * 3) The indices of the neighboring triangles are distinct for * each triangle. * 4) The indices of the edges are distinct in each triangle. * 5) Edge indices in ltri are in the range 0 to ne-1. * 6) Triangle indices in ltri, lvtx, and ledg are in the range * 0 to nt-1. * 7) Triangle ledg[ke] contains edge ke, and ledg[ke] is the * larger of the two relevant triangle indices for all ke. * 8) Triangle lvtx[kv] contains vertex kv for all kv. * 9) Adjacency information is consistent. * 10) Vertex coordinates vtxy are numbers (as opposed to Nan's). * 11) Triangle vertices are counterclockwise-ordered. * * On input: * * nc,nb,nv,nn,nt,ne = Numbers of boundary curves, boundary * edges, vertices, non-removable vertices, * triangles, and edges, respectively. * * ltri,lvtx,ledg,lnrv,lcrv,listn = Triangle list, vertex list, * edge list, non-removable vertex list, boundary curve * list, and list of non-removable vertex flags. * * If an integrity check is not satisfied, an error message is * printed, and the program is terminated. Otherwise a message * indicating success is printed. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void trmtst(int nc, int nb, int nv, int nn, int nt, int ne, double vtxy[][2], int ltri[][9], int lvtx[], int ledg[], int lnrv[], int lcrv[], unsigned char listn[]) { double sa, samin; int i, j, ke, kt, kt0, kv, kv0, kv1, kv2; unsigned char found; /* * Test for consistency of counts. */ if (nb < 3*nc || nn < 0 || nn > nv || ne != 3*nv-nb+3*nc-6 || nt != 2*nv-nb+2*nc-4) { printf("glmesh2 (trmtst): Inconsistent counts: nv = %d, " "nb = %d, nc = %d, nn = %d, nt = %d, ne = %d\n", nv, nb, nc, nn, nt, ne); exit(1); } /* * Test ltri values. */ for (kt = 0; kt < nt; kt++) { if (ltri[kt][0] < 0 || ltri[kt][0] >= nv || ltri[kt][1] < 0 || ltri[kt][1] >= nv || ltri[kt][2] < 0 || ltri[kt][2] >= nv || ltri[kt][3] < -1 || ltri[kt][3] >= nt || ltri[kt][4] < -1 || ltri[kt][4] >= nt || ltri[kt][5] < -1 || ltri[kt][5] >= nt || ltri[kt][6] < 0 || ltri[kt][6] >= ne || ltri[kt][7] < 0 || ltri[kt][7] >= ne || ltri[kt][8] < 0 || ltri[kt][8] >= ne) { printf("glmesh2 (trmtst): Invalid ltri[kt] entry: kt = " "%d\n", kt); exit(1); } if ( (ltri[kt][3] >= 0 && (ltri[kt][3] == ltri[kt][4] || ltri[kt][3] == ltri[kt][5])) || (ltri[kt][4] >= 0 && ltri[kt][4] == ltri[kt][5]) ) { printf("glmesh2 (trmtst): ltri[%d] has duplicate " "neighboring triangle indices\n", kt); exit(1); } if (ltri[kt][6] == ltri[kt][7] || ltri[kt][6] == ltri[kt][8] || ltri[kt][7] == ltri[kt][8]) { printf("glmesh2 (trmtst): ltri[%d] has duplicate " "edge indices\n", kt); exit(1); } } /* * Test ledg values. */ for (ke = 0; ke < ne; ke++) { kt = ledg[ke]; if (kt < 0 || kt >= nt) { printf("glmesh2 (trmtst): Invalid kt = ledg[ke] entry: " "ke = %d, kt = %d\n", ke, kt); exit(1); } if (ltri[kt][6] == ke) i = 3; else if (ltri[kt][7] == ke) i = 4; else if (ltri[kt][8] == ke) i = 5; else { printf("glmesh2 (trmtst): Invalid kt = ledg[ke] or " "ltri[kt] entry: ke = %d, kt = %d\n", ke, kt); exit(1); } kt0 = kt; kt = ltri[kt][i]; if (kt >= kt0) { printf("glmesh2 (trmtst): Invalid kt = ledg[ke] entry (not" " the smaller triangle index: ke = %d, kt = %d\n", ke, kt0); exit(1); } if (kt < 0) continue; if (ltri[kt][6] == ke) i = 3; else if (ltri[kt][7] == ke) i = 4; else if (ltri[kt][8] == ke) i = 5; else { printf("glmesh2 (trmtst): Invalid ledg[ke] or ltri[kt] " "entry: ke = %d, kt = %d\n", ke, kt); exit(1); } if (ltri[kt][i] != kt0) { printf("glmesh2 (trmtst): Inconsistency between ltri[kt1] " "and ltri[kt2], kt1 = %d, kt2 = %d\n", kt0, kt); exit(1); } } /* * Test lvtx values. */ for (kv = 0; kv < nv; kv++) { kt = lvtx[kv]; if (kt < 0 || kt >= nt) { printf("glmesh2 (trmtst): Invalid kt = lvtx[kv] entry: " "kv = %d, kt = %d\n", kv, kt); exit(1); } if (ltri[kt][0] == kv) i = 0; else if (ltri[kt][1] == kv) i = 1; else if (ltri[kt][2] == kv) i = 2; else { printf("glmesh2 (trmtst): Invalid kt = lvtx[kv] or " "ltri[kt] entry: kv = %d, kt = %d\n", kv, kt); exit(1); } } /* * Test lnrv values and flags in listn. */ for (i = 0; i < nn; i++) { kv = lnrv[i]; if (kv < 0 || kv >= nv) { printf("glmesh2 (trmtst): Invalid kv = lnrv[k] entry: k =" " %d, kv = %d\n", i, kv); exit(1); } kt = lvtx[kv]; if (ltri[kt][0] == kv) j = 2; else if (ltri[kt][1] == kv) j = 0; else j = 1; if (ltri[kt][j+3] >= 0) { printf("glmesh2 (trmtst): Vertex kv = lnrv[k] is not a " "boundary vertex: k = %d, kv = %d\n", i, kv); exit(1); } if (!listn[kv]) { printf("glmesh2 (trmtst): Vertex %d is in lnrv but has " "zero flag value in listn\n", kv); exit(1); } } for (kv = 0; kv < nv; kv++) { if (!listn[kv]) continue; found = 0; for (i = 0; i < nn; i++) { if (lnrv[i] == kv) { found = 1; break; } } if (!found) { printf("glmesh2 (trmtst): Vertex %d has listn flag value " "1, but is not in lnrv\n", kv); exit(1); } } /* * Test lcrv values. */ for (i = 0; i < nc; i++) { kv = lcrv[i]; if (kv < 0 || kv >= nv) { printf("glmesh2 (trmtst): Invalid kv = lcrv[k] entry: k =" " %d, kv = %d\n", i, kv); exit(1); } kt = lvtx[kv]; if (ltri[kt][0] == kv) j = 2; else if (ltri[kt][1] == kv) j = 0; else j = 1; if (ltri[kt][j+3] >= 0) { printf("glmesh2 (trmtst): Vertex kv = lcrv[k] is not a " "boundary vertex: k = %d, kv = %d\n", i, kv); exit(1); } } /* * Test vtxy values. */ for (i = 0; i < nv; i++) { if (vtxy[i][0] != vtxy[i][0] || vtxy[i][1] != vtxy[i][1]) { printf("glmesh2 (trmtst): Vertex kv has an invalid " "component (Nan): kv = %d\n", i); exit(1); } } /* * Compute the minimum signed triangle area. A negative value * implies a triangle with clockwise-ordered vertices. */ samin = 1.0e37; for (kt = 0; kt < nt; kt++) { kv0 = ltri[kt][0]; kv1 = ltri[kt][1]; kv2 = ltri[kt][2]; sa = (vtxy[kv1][0]-vtxy[kv0][0])*(vtxy[kv2][1]-vtxy[kv0][1]) - (vtxy[kv2][0]-vtxy[kv0][0])*(vtxy[kv1][1]-vtxy[kv0][1]); if (sa < samin) samin = sa; } samin *= 0.5; printf("glmesh2 (trmtst): Minimum signed triangle area = %e\n", samin); if (samin < 0.0) return; /* * No error encountered. */ printf("Triangulation integrity tests passed.\n"); return; } /* End of trmtst */ /* *-------------------------------------------------------------------- * * trprint: Print the triangle mesh data structure, along with * triangle shape-quality statistics, on the standard * output unit. * * On input: * * nc,nb,nv,nn,nt,ne = Numbers of boundary curves, boundary * edges, vertices, non-removable vertices, * triangles, and edges, respectively. * * vtxy = Type double array dimensioned nv by 2 containing the * Cartesian coordinates of the vertices. * * uv = Type double array dimensioned nv by 2 containing pairs * of function values (u,v) at the vertices. * * ltri,lvtx,ledg,lnrv,lcrv = Triangle list, vertex list, edge * list, non-removable vertex list, * and boundary curve list. * * glmesh2 functions called: initR, trqual1, trqual2 * *-------------------------------------------------------------------- */ static void trprint(int nc, int nb, int nv, int nn, int nt, int ne, double vtxy[][2], double uv[][2], int ltri[][9], int lvtx[], int ledg[], int lnrv[], int lcrv[]) { char *s; char *separator = "\n________________________________________" "________________________________________\n"; double cfac, emin, pfair, pgood, q, qmax, qmean, qmin, scamx; int i, imin, kt, kt0, kv, kv1, kv2, kv3, nnb; /* * Store listR if necessary. */ if (newR) initR(); /* * Print title and counts. */ printf("%s", separator); printf(" Triangle Mesh Data Structure\n\n" " Nc = %d boundary curves\n" " Nb = %d boundary edges (boundary vertices)\n" " Nv = %d vertices\n" " Nn = %d non-removable boundary vertices\n" " Nt = %d triangles\n" " Ne = %d edges\n\n", nc, nb, nv, nn, nt, ne); /* * Print the vertex coordinates, (u,v) pairs, vertex list lvtx, * and *'s after the non-removable vertices. */ printf(" Vertex List\n\n" " Vertex x y u v" " lvtx Non-removable\n\n"); for (kv = 0; kv < nv; kv++) { s = ""; for (i = 0; i < nn; i++) { if (lnrv[i] == kv) { s = "*"; break; } } printf(" %6.1d %10.3e %10.3e %10.3e %10.3e %6.1d" " %s\n", kv, vtxy[kv][0], vtxy[kv][1], uv[kv][0], uv[kv][1], lvtx[kv], s); } /* * Print the triangle list. */ printf("\n\n Triangle List\n\n" "Triangle Vertex Indices Neighboring Triangles " "Edge Indices\n\n"); for (i = 0; i < nt; i++) printf(" %6.1d %6.1d %6.1d %6.1d %6.1d %6.1d %6.1d" " %6.1d %6.1d %6.1d\n", i, ltri[i][0], ltri[i][1], ltri[i][2], ltri[i][3], ltri[i][4], ltri[i][5], ltri[i][6], ltri[i][7], ltri[i][8]); /* * Print the edge list. */ printf("\n\n Edge list (triangle index for each edge):\n\n"); for (i = 0; i < ne-8; i += 9) printf(" %6.1d %6.1d %6.1d %6.1d %6.1d %6.1d %6.1d" " %6.1d %6.1d\n", ledg[i], ledg[i+1], ledg[i+2], ledg[i+3], ledg[i+4], ledg[i+5], ledg[i+6], ledg[i+7], ledg[i+8]); for (i = ne-(ne%9); i < ne; i++) printf(" %6.1d", ledg[i]); printf("\n"); /* * Print the boundary curve list. */ printf("\n\n Boundary curve list (boundary vertex index for each" " curve):\n\n"); for (i = 0; i < nc-8; i += 9) printf(" %6.1d %6.1d %6.1d %6.1d %6.1d %6.1d %6.1d" " %6.1d %6.1d\n", lcrv[i], lcrv[i+1], lcrv[i+2], lcrv[i+3], lcrv[i+4], lcrv[i+5], lcrv[i+6], lcrv[i+7], lcrv[i+8]); for (i = nc-(nc%9); i < nc; i++) printf(" %6.1d", lcrv[i]); printf("\n"); /* * Compute and print shape-quality statistics. * * q1: Area/(Mean squared side length). */ qmin = 2.0; qmean = 0.0; qmax = 0.0; pfair = 0.0; pgood = 0.0; for (kt = 0; kt < nt; kt++) { kv1 = ltri[kt][0]; kv2 = ltri[kt][1]; kv3 = ltri[kt][2]; q = trqual1(vtxy[kv1][0], vtxy[kv1][1], vtxy[kv2][0], vtxy[kv2][1], vtxy[kv3][0], vtxy[kv3][1]); if (q < qmin) qmin = q; qmean += q; if (q > qmax) qmax = q; if (q > 0.6) pfair += 1.0; if (q > 0.866) pgood += 1.0; } qmean /= (GLdouble) nt; pfair = 100.0*pfair/(GLdouble) nt; pgood = 100.0*pgood/(GLdouble) nt; printf("\n Triangle shape quality q1 = " "Area/Mean-squared-side-length (normalized):\n"); printf(" q1_min = %5.3f Percent poor (q1 < 0.6) = %6.2f\n", qmin, 100.0-pfair); printf(" q1_mean = %5.3f Percent fair (q1 >= 0.6) = %6.2f\n", qmean, pfair); printf(" q1_max = %5.3f Percent good (q1 >= 0.866) = %6.2f\n", qmax, pgood); /* * q2: Inradius/Circumradius. */ qmin = 2.0; qmean = 0.0; qmax = 0.0; pfair = 0.0; pgood = 0.0; for (kt = 0; kt < nt; kt++) { kv1 = ltri[kt][0]; kv2 = ltri[kt][1]; kv3 = ltri[kt][2]; q = trqual2(vtxy[kv1][0], vtxy[kv1][1], vtxy[kv2][0], vtxy[kv2][1], vtxy[kv3][0], vtxy[kv3][1]); if (q < qmin) qmin = q; qmean += q; if (q > qmax) qmax = q; if (q > 0.5) pfair += 1.0; if (q > 0.7) pgood += 1.0; } qmean /= (GLdouble) nt; pfair = 100.0*pfair/(GLdouble) nt; pgood = 100.0*pgood/(GLdouble) nt; printf("\n Triangle shape quality q2 = " "Inradius/Circumradius (normalized):\n"); printf(" q2_min = %5.3f Percent poor (q2 < 0.5) = %6.2f\n", qmin, 100.0-pfair); printf(" q2_mean = %5.3f Percent fair (q2 >= 0.5) = %6.2f\n", qmean, pfair); printf(" q2_max = %5.3f Percent good (q2 >= 0.7) = %6.2f\n", qmax, pgood); /* * Maximum and mean absolute deviation of vertex degree from 6: */ qmax = 0.0; qmean = 0.0; for (kv = 0; kv < nv; kv++) { kt0 = lvtx[kv]; kt = kt0; nnb = 0; do { if (ltri[kt][0] == kv) i = 1; else if (ltri[kt][1] == kv) i = 2; else i = 0; kt = ltri[kt][i+3]; nnb++; } while (kt >= 0 && kt != kt0); if (kt < 0) continue; q = (double) abs(nnb-6); qmean += q; if (q > qmax) qmax = q; } qmean /= (GLdouble) (nv-nb); printf("\n Topological triangle shape quality " "(degrees of interior vertices):\n"); printf(" mean absolute deviation from 6 = %5.1f\n", qmean); printf(" maximum absolute deviation from 6 = %5.1f\n", qmax); /* * Minimum angle. */ qmin = 180.0; qmean = 0.0; qmax = 0.0; pgood = 0.0; cfac = 180.0/acos(-1.0); /* Radians-to-degrees scale factor */ for (kt = 0; kt < nt; kt++) { badTriangle(kt, &imin, &emin, &scamx); q = cfac*acos(sqrt(scamx)); if (q < qmin) qmin = q; qmean += q; if (q > qmax) qmax = q; if (q > min_angle) pgood += 1.0; } qmean /= (GLdouble) nt; pgood = 100.0*pgood/(GLdouble) nt; printf("\n Minimum angle a:\n"); printf(" a_min = %5.1f\n", qmin); printf(" a_mean = %5.1f\n", qmean); printf(" a_max = %5.1f\n", qmax); printf(" Percent good (a >= %d) = %6.2f\n", min_angle, pgood); printf("%s", separator); fflush(stdout); return; } /* End of trprint */ /* *-------------------------------------------------------------------- * * trqual1: Compute a measure of triangle quality, * * q = 4*sqrt(3)*A/(s1*s1 + s2*s2 + s3*s3), * * where A is the triangle area, and s1, s2, and s3 are * the side lengths. The scaling is chosen so that an * equilateral triangle has q = 1, which is optimal. * * We define three categories: * * 1) Good: q >= sqrt(3)/2. The largest angle is at most pi/2, and * the smallest angle is at least arccos(.8) in this case. * * 2) Fair: 0.6 <= q < sqrt(3)/2. The largest angle is at most * 2*pi/3, and the smallest angle is at least * arccos(13/14) in this case. * * 3) Poor: q < 0.6. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static double trqual1(double x1, double y1, double x2, double y2, double x3, double y3) { double a, dx1, dx2, dx3, dy1, dy2, dy3, s1, s2, s3; /* * Compute squared side lengths s1, s2, and s3. */ dx1 = x3-x2; dy1 = y3-y2; dx2 = x1-x3; dy2 = y1-y3; dx3 = x2-x1; dy3 = y2-y1; s1 = dx1*dx1 + dy1*dy1; s2 = dx2*dx2 + dy2*dy2; s3 = dx3*dx3 + dy3*dy3; /* * Compute a = 2*area. */ a = fabs( dx3*dy2 - dx2*dy3 ); return (2.0*sqrt(3.0)*a/(s1 + s2 + s3)); } /* End of trqual1 */ /* *-------------------------------------------------------------------- * * trqual2: Compute a measure of triangle quality, * * q = 2*r/R = (s2+s3-s1)*(s3+s1-s2)*(s1+s2-s3)/s1*s2*s3, * * where r and R denote the inradius and circumradius of * the triangle, and s1, s2, and s3 are the side lengths. * The value ranges from 0 for a degenerate triangle to 1 * for an equilateral triangle. * * We define three categories: * * 1) Good: q >= 0.7 * * 2) Fair: 0.5 <= q < 0.7. * * 3) Poor: q < 0.5. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static double trqual2(double x1, double y1, double x2, double y2, double x3, double y3) { double dx1, dx2, dx3, dy1, dy2, dy3, s1, s2, s3; /* * Compute side lengths s1, s2, and s3. */ dx1 = x3-x2; dy1 = y3-y2; dx2 = x1-x3; dy2 = y1-y3; dx3 = x2-x1; dy3 = y2-y1; s1 = sqrt(dx1*dx1 + dy1*dy1); s2 = sqrt(dx2*dx2 + dy2*dy2); s3 = sqrt(dx3*dx3 + dy3*dy3); return ((s2+s3-s1)*(s3+s1-s2)*(s1+s2-s3)/(s1*s2*s3)); } /* End of trqual2 */ /* *-------------------------------------------------------------------- * * trwrite: Output the triangle mesh to a file. * * On input: * * fname = File name for output. * * nc,nb,nv,nt,ne = Numbers of boundary curves, boundary edges, * vertices, triangles, and edges, respectively. * * vtxy = Type double array dimensioned nv by 2 containing the * Cartesian coordinates of the vertices. * * uv = Type double array dimensioned nv by 2 containing pairs * of function values (u,v) at the vertices. * * ltri = Triangle list. * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void trwrite(char *fname, int nc, int nb, int nv, int nt, int ne, double vtxy[][2], double uv[][2], int ltri[][9]) { FILE *fpw; /* File pointer */ time_t curtime; /* Current calendar time */ struct tm *loctime; /* Local time representation */ int linsiz = 160; /* Buffer size */ char line[160]; /* Buffer for output line */ int kt, kv, nbc = 0; /* * Open file fname with write access. */ fpw = fopen(fname, "w"); if (fpw == NULL) { printf("glmesh2: Cannot open file %s.\n", fname); exit(1); } /* * Output header comments. */ curtime = time (NULL); loctime = localtime (&curtime); strftime(line,linsiz,"#\n# Data file created by glmesh2 on " "%b %d, %Y\n", loctime); fputs(line,fpw); fprintf(fpw,"# Nc = %d, Nb = %d, Nv = %d, Nt = %d, Ne = %d\n", nc, nb, nv, nt, ne); fprintf(fpw,"#\n# Vertex list: Nv followed by Nv (x,y,u,v) " "values\n"); /* * Output vertex list. */ fprintf(fpw," %d\n", nv); for (kv = 0; kv < nv; kv++) { fprintf(fpw," %19.12e %19.12e %19.12e %19.12e\n", vtxy[kv][0], vtxy[kv][1], uv[kv][0], uv[kv][1]); } /* * Output the triangle list. */ fprintf(fpw,"#\n# Triangle list: Nt followed by Nt CCW-ordered" " triples of vertex indices\n"); fprintf(fpw," %d\n", nt); for (kt = 0; kt < nt; kt++) { fprintf(fpw," %9d %9d %9d\n", ltri[kt][0], ltri[kt][1], ltri[kt][2]); } /* * Output a boundary segment list. */ fprintf(fpw,"#\n# Boundary segment list: Nbc followed by Nbc" " triples (type,ke1,ke2)\n"); fprintf(fpw," %d\n", nbc); return; } /* End of trwrite */ /* *-------------------------------------------------------------------- * * unrefine: Unrefine the mesh by reversing the previous refinement * by refineRd, refineRe, or refineRt unless there was no * previous refinement or the previous refinement was * rendered non-reversible by deletion of a vertex (in * Function coarsenRv or removeRt). * * Note that the vertices inserted by the previous refinement will * not all be removed if some were converted to non-removable (by * Function convertRv) since the refinement. Unrefinement does not * restore the previous mesh in this case. Also, the Delaunay edge * swaps associated with refineRd are not reversed. Furthermore, if * an unrefinement fails to remove all vertices, a subsequent * unrefinement will not reverse the original refinement. * * glmesh2 function called: delvtx * *-------------------------------------------------------------------- */ static void unrefine(void) { int kv, nfail, nv0; if (nref <= 0) { printf("\n Previous refinement cannot be reversed.\n"); return; } /* * Pop rstack into nv0. */ nref--; nv0 = rstack[nref]; /* * Delete vertex sequence nv-1, nv-2, ..., nv0. */ nfail = 0; for (kv = nv-1; kv >= nv0; kv--) { if (!delvtx(kv)) nfail++; } if (nfail) { printf("\n Unrefinement failed to remove all vertices.\n"); } return; } /* End of unrefine */ /* *-------------------------------------------------------------------- * * unswap: Reverse the effect of a call to Function swap: swap * diagonal edges in a convex quadrilateral defined by a * pair of adjacent triangles, with the choice of new * triangle indices reversed from that of Function swap, * so that the application of either function followed by * the other leaves the triangle mesh unaltered, not only * in terms of geometry and topology, but also in the data * structure. Applying the same function twice has the * effect of reversing the indices of the two triangles, * and altering the cyclical ordering of vertices. That has * side effects on other data structures that use triangle * indices. * * On input: * * ke = Index of the edge to be swapped. * * It is assumed that the quadrilateral is convex. No variables are * altered if edge ke is not shared by two triangles. * * Global variables required: ledg, ltri, lvtx * * glmesh2 functions called: None * *-------------------------------------------------------------------- */ static void unswap(int ke) { int i, j, ken1, ken2, kt1, kt2, ktn1, ktn2, kv1, kv2, kvn1, kvn2; /* * Set kt1 and kt2 to the adjacent triangles with shared edge ke, * set (kv1,kvn1,kvn2) and (kv2,kvn2,kvn1) to the vertex indices of * triangles kt1 and kt2, respectively, and set i and j to the column * indices of kv1 in kt1 and kv2 in kt2, respectively. */ kt1 = ledg[ke]; if (ltri[kt1][6] == ke) i = 0; else if (ltri[kt1][7] == ke) i = 1; else i = 2; kt2 = ltri[kt1][i+3]; if (kt2 < 0) return; if (ltri[kt2][6] == ke) j = 0; else if (ltri[kt2][7] == ke) j = 1; else j = 2; kv1 = ltri[kt1][i]; kv2 = ltri[kt2][j]; kvn1 = ltri[kt1][(i+1)%3]; kvn2 = ltri[kt2][(j+1)%3]; ktn1 = ltri[kt2][(j+2)%3+3]; ktn2 = ltri[kt1][(i+2)%3+3]; ken1 = ltri[kt2][(j+2)%3+6]; ken2 = ltri[kt1][(i+2)%3+6]; /* * If i = j = 0, the initial ltri rows are * * kt1 = (kv1,kvn1,kvn2,kt2,ktn2,x,ke,ken2,x) * kt2 = (kv2,kvn2,kvn1,kt1,ktn1,x,ke,ken1,x) * * and the new rows, following the swap are * * kt1 = (kv1,kv2,kvn2,ktn1,x,kt2,ken1,x,ke) * kt2 = (kv2,kv1,kvn1,ktn2,x,kt1,ken2,x,ke) * * where the entries labeled x retain their values and need not be * referenced. */ ltri[kt1][(i+1)%3] = kv2; ltri[kt2][(j+1)%3] = kv1; ltri[kt1][i+3] = ktn1; ltri[kt2][j+3] = ktn2; ltri[kt1][(i+2)%3+3] = kt2; ltri[kt2][(j+2)%3+3] = kt1; ltri[kt1][i+6] = ken1; ltri[kt2][j+6] = ken2; ltri[kt1][(i+2)%3+6] = ke; ltri[kt2][(j+2)%3+6] = ke; /* * Correct neighboring triangle ltri entries if necessary. */ if (ktn1 >= 0) { if (ltri[ktn1][3] == kt2) ltri[ktn1][3] = kt1; else if (ltri[ktn1][4] == kt2) ltri[ktn1][4] = kt1; else ltri[ktn1][5] = kt1; } if (ktn2 >= 0) { if (ltri[ktn2][3] == kt1) ltri[ktn2][3] = kt2; else if (ltri[ktn2][4] == kt1) ltri[ktn2][4] = kt2; else ltri[ktn2][5] = kt2; } /* * Adjust lvtx entries. */ if (lvtx[kv1] == kt1) lvtx[kv1] = kt2; if (lvtx[kv2] == kt2) lvtx[kv2] = kt1; if (lvtx[kvn1] == kt1) lvtx[kvn1] = kt2; if (lvtx[kvn2] == kt2) lvtx[kvn2] = kt1; /* * Adjust ledg entries, preserving the property that the larger of * the two triangle indices is used. */ ledg[ken1] = kt1; if (ktn1 > kt1) ledg[ken1] = ktn1; ledg[ken2] = kt2; if (ktn2 > kt2) ledg[ken2] = ktn2; return; } /* End of unswap */ /* *-------------------------------------------------------------------- * * update: Update a Delaunay triangulation in the vicinity of a * newly inserted vertex v by applying a sequence of swap * tests and the appropriate swaps to the edges opposite v. * If an edge is swapped out, there are two new edges that * must be tested for swaps. This is the essence of * Lawson's incremental algorithm for constructing a * Delaunay triangulation. Also, the boundary edges * opposite v, if contained in R, are tested for encroach- * ment by v, and added to Qe if encroached. Optionally, * the bad triangles containing v, if in R, are added to Qt. * * On input: kv is the index (0 to nv-1) of the vertex v whose star * (set of containing triangles) is to be updated. * * tritest is a flag with value 1 iff the triangles that * contain v are to be tested by Function badTriangle for * inclusion in Qt. * * On output: ns is the number of swaps, and inde contains the * indices of the swapped edges, in the order in which * the swaps occurred, in the first ns positions. The * input triangulation can be restored by swapping edges * in the reverse order: unswap(inde[k]) for k = ns-1, * ns-2, ..., 1, 0. * * Reference: C. L. Lawson (1977), Software for C^1 surface interp- * olation, in Mathematical Software III (J. R. Rice, * ed.), Academic Press, New York, pp. 161-194. * * glmesh2 functions called: badTriangle, encroached, insertQt, * swap, swptst * *-------------------------------------------------------------------- */ static void update(int kv, unsigned char tritest, int *ns, int inde[]) { double emin, scamx; int i, imin, in1, in2, io1, io2, j, ke, kt, kt0, ktn, kv1, kv2, ls; unsigned char enq, swp; /* * Initialization: ls = length of inde. */ ls = 0; kt0 = lvtx[kv]; kt = kt0; /* * Loop on triangles kt containing kv with opposite edge ke. * enq is a flag with value 1 iff ke is added to Qe. */ do { enq = 0; if (ltri[kt][0] == kv) i = 0; else if (ltri[kt][1] == kv) i = 1; else i = 2; in1 = ltri[kt][i]; io1 = ltri[kt][(i+1)%3]; io2 = ltri[kt][(i+2)%3]; ktn = ltri[kt][i+3]; ke = ltri[kt][i+6]; if (ktn >= 0) { if (ltri[ktn][6] == ke) j = 0; else if (ltri[ktn][7] == ke) j = 1; else j = 2; in2 = ltri[ktn][j]; } else if (listR[io1] && listR[io2]) { /* * Test boundary edge ke for encroachment. kv1 and kv2 are unused * copies of io1 and io2. */ if (encroached(ke, &kv1, &kv2)) { /* * Test for sufficient storage in Qe, and enqueue ke. */ if (neq == neqmax) reallocQe(); enq = 1; qeback = (qeback+1)%neqmax; Qe[qeback].index_e = ke; Qe[qeback].index_v1 = io1; Qe[qeback].index_v2 = io2; neq++; } } swp = (ktn >= 0 && swptst(in1,in2,io1,io2)); if (swp) { swap(ke); inde[ls] = ke; ls++; } else { if (tritest && !enq && listR[kv] && listR[io1] && listR[io2]) { /* * Add triangle kt to the priority queue Qt if it is bad. */ if (badTriangle(kt, &imin, &emin, &scamx)) { insertQt(kt, imin, emin); } } /* * Advance to the next triangle. */ i = (i+1)%3; kt = ltri[kt][i+3]; } } while (kt >= 0 && (kt != kt0 || swp)); *ns = ls; return; } /* End of update */ /* *-------------------------------------------------------------------- * * xlate: Translate glut window coordinates (wx,wy), with depth wz, * to object coordinates (x,y,z) by inverting the mappings * associated with the modelview, projection, and viewport * matrices. * * On input: wx and wy are the integer window-relative coordinates * in the range 0 to w-1 and 0 to h-1, respectively, for * window width w and height h, and wz is the depth in * the range 0 to 1.0. The origin is at the upper left * corner of the window. * * On output: x, y, and z are the objects coordinates. * * Global variables required: None * * glmesh2 functions required: None * *-------------------------------------------------------------------- */ static void xlate(GLint wx, GLint wy, GLdouble wz, GLdouble *x, GLdouble *y, GLdouble *z) { GLdouble mvmatrix[16]; /* Modelview matrix */ GLdouble projmatrix[16]; /* Projection matrix */ GLint viewport[4]; /* Components of Viewport mapping */ GLint yr; /* OpenGL y coordinate */ glGetIntegerv(GL_VIEWPORT, viewport); glGetDoublev(GL_MODELVIEW_MATRIX, mvmatrix); glGetDoublev(GL_PROJECTION_MATRIX, projmatrix); /* * Compute yr = h-1-wy for window height h. */ yr = glutGet(GLUT_WINDOW_HEIGHT) - 1 - wy; gluUnProject((GLdouble) wx, (GLdouble) yr, wz, mvmatrix, projmatrix, viewport, x, y, z); return; } /* End of xlate */ /* *-------------------------------------------------------------------- * * main: Input data, create windows, create a menu, register call- * back functions, and enter a loop. * * glmesh2 functions called: getCorners, init0, inputData, makeMenu, * storeV * * glmesh2 callback functions registered: display, dlgDisplay, * dlgkey, dlgMouse, dlgReshape, key, motion, mouse, * passiveMotion, reshape, specialKey * *-------------------------------------------------------------------- */ int main(int argc, char *argv[]) { inputData(argc, argv); /* Input data */ storeV(); /* Initialize global parameters */ getCorners(); /* Initialize lnrv to corners */ glutInit(&argc, argv); /* Initialize Glut library */ /* * Create the primary window (surface display): window 0. */ glutInitWindowPosition(50, 50); glutInitWindowSize(500, 500); /* Size in pixels */ glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE); window[0] = glutCreateWindow(argv[1]); /* Create window 0 */ makeMenu(); /* Create a menu attached to the right button */ glutKeyboardFunc(key); /* Register callbacks */ glutSpecialFunc(specialKey); glutMouseFunc(mouse); glutMotionFunc(motion); glutPassiveMotionFunc(passiveMotion); glutReshapeFunc(reshape); glutDisplayFunc(display); init0(); /* One-time initialization */ /* * Create the secondary window (dialog): window 1. */ glutInitWindowPosition(500, 250); glutInitWindowSize(400, 200); /* Size in pixels */ glutInitDisplayMode(GLUT_RGB | GLUT_SINGLE); window[1] = glutCreateWindow("Dialog"); /* Create window 1 */ glClearColor (0.5, 0.5, 0.5, 1.0); /* Grey background */ glutKeyboardFunc(dlgKey); /* Register callbacks */ glutMouseFunc(dlgMouse); glutReshapeFunc(dlgReshape); glutDisplayFunc(dlgDisplay); glutHideWindow(); /* Initially hidden */ glutSetWindow(window[0]); /* Make window 0 current */ glutMainLoop(); return 0; } /* End of main */