Adding the tetra-metrics functionality, we compute the quality and ed…

src/geometry.c

@@ -1,6 +1,7 @@
 #include "medit.h"
 #include "extern.h"
 #include "sproto.h"
+#include "mesh.h"
 
 extern int getIso(pScene sc,double norm,double *hsv);
 

src/grafic.h

@@ -1,6 +1,8 @@
 #ifndef _GRAFIC_H
 #define _GRAFIC_H
 
+#include "mesh.h"
+
 #define MAX_LIST  4
 #define MAXISO    5
 
@@ -232,6 +234,12 @@ typedef struct scene {
   ubyte      type;
   ubyte      isotyp;
   ubyte      picked;
+
+  /* Adding the field for the quality mode of the tetrahedra */
+  TetraMetrics currentTetraMetrics;  // Store the metrics of the selected tetrahedron
+  ubyte displayTetraMetrics;         // Flag activating/deactivating the 
+				     // display of the tetrahedron metrics
+
 } Scene;
 typedef Scene * pScene;
 

src/medit.h

@@ -88,4 +88,4 @@ typedef Canvas * pCanvas;
 
 
 
-#endif
+#endif

src/mesh.c

@@ -1,6 +1,7 @@
 #include "medit.h"
 #include "extern.h"
 #include "sproto.h"
+#include "mesh.h" 
 
 #define FLOAT_MAX  1.e20
 ubyte   dosurf;
@@ -475,3 +476,92 @@ int meshUpdate(pScene sc,pMesh mesh) {
 
   return(1);
 }
+
+// tolerance to prevent a division by zero
+static const double epsd2 = 1e-12;
+
+TetraMetrics compute_mesh_quality_edges_length(pMesh mesh, int tetraIndex) {
+    TetraMetrics metrics;
+
+    // Accessing the targeted tetrahedron
+    pTetra tet = &mesh->tetra[tetraIndex];
+
+    // Getting the 4 points forming the tetrahedron
+    pPoint p1 = &mesh->point[tet->v[0]];
+    pPoint p2 = &mesh->point[tet->v[1]];
+    pPoint p3 = &mesh->point[tet->v[2]];
+    pPoint p4 = &mesh->point[tet->v[3]];
+
+    // Computing the 6 edge lengths
+    double d12 = sqrt(pow(p2->c[0] - p1->c[0], 2) +
+                      pow(p2->c[1] - p1->c[1], 2) +
+                      pow(p2->c[2] - p1->c[2], 2));
+    double d13 = sqrt(pow(p3->c[0] - p1->c[0], 2) +
+                      pow(p3->c[1] - p1->c[1], 2) +
+                      pow(p3->c[2] - p1->c[2], 2));
+    double d14 = sqrt(pow(p4->c[0] - p1->c[0], 2) +
+                      pow(p4->c[1] - p1->c[1], 2) +
+                      pow(p4->c[2] - p1->c[2], 2));
+    double d23 = sqrt(pow(p3->c[0] - p2->c[0], 2) +
+                      pow(p3->c[1] - p2->c[1], 2) +
+                      pow(p3->c[2] - p2->c[2], 2));
+    double d24 = sqrt(pow(p4->c[0] - p2->c[0], 2) +
+                      pow(p4->c[1] - p2->c[1], 2) +
+                      pow(p4->c[2] - p2->c[2], 2));
+    double d34 = sqrt(pow(p4->c[0] - p3->c[0], 2) +
+                      pow(p4->c[1] - p3->c[1], 2) +
+                      pow(p4->c[2] - p3->c[2], 2));
+
+    // Storing the 6 edge lengths in our structure
+    metrics.edgeLengths[0] = d12;
+    metrics.edgeLengths[1] = d13;
+    metrics.edgeLengths[2] = d14;
+    metrics.edgeLengths[3] = d23;
+    metrics.edgeLengths[4] = d24;
+    metrics.edgeLengths[5] = d34;
+
+    // Computing the volume using 3 arrays coming from the point p1 :
+    // ab = p2 - p1, ac = p3 - p1, ad = p4 - p1
+    double ab[3] = { p2->c[0] - p1->c[0],
+                     p2->c[1] - p1->c[1],
+                     p2->c[2] - p1->c[2] };
+    double ac[3] = { p3->c[0] - p1->c[0],
+                     p3->c[1] - p1->c[1],
+                     p3->c[2] - p1->c[2] };
+    double ad[3] = { p4->c[0] - p1->c[0],
+                     p4->c[1] - p1->c[1],
+                     p4->c[2] - p1->c[2] };
+
+    // Computing the cross product (ac x ad)
+    double cross[3] = {
+        ac[1]*ad[2] - ac[2]*ad[1],
+        ac[2]*ad[0] - ac[0]*ad[2],
+        ac[0]*ad[1] - ac[1]*ad[0]
+    };
+
+    // Computing the scalar product (ab · (ac x ad)) to get the 6*Volume
+    double vol6 = ab[0]*cross[0] + ab[1]*cross[1] + ab[2]*cross[2];
+
+    // Verifying the volume is not too small 
+    // (to prevent the division by zero)
+    if(fabs(vol6) < epsd2) {
+        metrics.quality = 0.0;
+        return metrics;
+    }
+
+    // Computing S = sum of the squared edge lengths
+    double S = d12*d12 + d13*d13 + d14*d14 + d23*d23 + d24*d24 + d34*d34;
+    if(S < epsd2) {
+        metrics.quality = 0.0;
+        return metrics;
+    }
+
+    // Computing the quality with respect to the modified formula :
+    // For a regular tetrahedron of edge length l, 
+    // The ratio computed is 1/(12sqrt(3)), we multiply 
+    // by 12sqrt(3) to obtain 1.
+    metrics.quality = (vol6 * (12 * sqrt(3))) / (S * sqrt(S));  
+
+    return metrics;
+}
+

src/mesh.h

@@ -16,6 +16,17 @@ typedef unsigned short uShort;
 #endif
 
 
+#ifndef _TETRA_METRICS_H
+#define _TETRA_METRICS_H
+
+typedef struct {
+    int id;                // Id of the tetrahedron	
+    double quality;        // Quality of the tetrahedron, with tet_quality
+    double edgeLengths[6]; // The lengths of the 6 tetrahedron edges
+} TetraMetrics;
+
+#endif
+
 typedef struct spoint {
   double    c[3];
   int       tmp,mark;

src/mouse.c

@@ -1,6 +1,7 @@
 #include "medit.h"
 #include "extern.h"
 #include "sproto.h"
+#include "mesh.h"
 
 #ifndef  ON
 #define  ON     1
@@ -150,20 +151,26 @@ void mouse(int button,int state,int x,int y) {
   tr->mstate  = state;
   tr->mbutton = button;
 
-  /* check if ctrl-shift-alt pressed */
+  /* check the actions on the keys */
   keyact = glutGetModifiers();
 
   if ( state == GLUT_DOWN ) {
     tracking = GL_TRUE;
     lasttime = glutGet(GLUT_ELAPSED_TIME);
 
+  /* check if shift+click pressed */
     if ( button == GLUT_LEFT_BUTTON ) {
       if ( keyact & GLUT_ACTIVE_SHIFT ) {
-        /* entity designation */
-        picking = GL_TRUE;
-				if ( sc->picklist ) glDeleteLists(sc->picklist,1);
-				sc->picklist = pickingScene(sc,x,y,0);
-				return;
+	      // Picking mode for computing tetrahedron metrics using SHIFT.
+	      picking = GL_TRUE;
+
+	      // Clean up any previous display list
+	      if ( sc->picklist ) glDeleteLists(sc->picklist,1);
+
+	      // Create a new display list to highlight the selected object
+	      sc->picklist = pickingScene(sc,x,y,0);
+	      return;
+
       }
       else if ( keyact & GLUT_ACTIVE_ALT ) {
 	    	/* zoom */

src/picking.c

@@ -1,7 +1,7 @@
 #include "medit.h"
 #include "extern.h"
 #include "sproto.h"
-
+#include "mesh.h"
 
 typedef struct color {
   GLuint  rMask,gMask,bMask,aMask;
@@ -899,6 +899,7 @@ GLuint pickingScene(pScene sc,int x,int y,int ident) {
     case LTets:
       drawTets(sc,mesh,refitem);
       infoEntity(sc,mesh,refitem,LTets);
+      pickingTetrahedron(sc, x, y);
       break;
     case LHexa:
       drawHexa(sc,mesh,refitem);
@@ -915,6 +916,38 @@ GLuint pickingScene(pScene sc,int x,int y,int ident) {
   return(dlist); 
 }
 
+int pickingTetrahedron(pScene sc, int x, int y) {	
+
+
+    // reftype/refitem have just been updated
+    if ( reftype == LTets ){
+	int tetIndex = refitem; 
+
+	// Computing and storing the metrics
+	TetraMetrics metrics = compute_mesh_quality_edges_length(cv.mesh[sc->idmesh], tetIndex);
+	metrics.id = tetIndex;
+	sc->currentTetraMetrics = metrics;
+	sc->displayTetraMetrics = 1;
+
+	// Print the computed metrics in the shell
+	printf("  Quality: %g\n", metrics.quality);
+	printf("  Edge Lengths: %.5g, %.5g, %.5g, %.5g, %.5g, %.5g\n",
+			metrics.edgeLengths[0],
+			metrics.edgeLengths[1],
+			metrics.edgeLengths[2],
+			metrics.edgeLengths[3],
+			metrics.edgeLengths[4],
+			metrics.edgeLengths[5]);
+
+	return tetIndex;
+
+    }      
+
+    // Not a tetrahedron
+    sc->displayTetraMetrics = 0;
+    return -1;
+
+}
 
 GLuint pickItem(pMesh mesh,pScene sc,int numit) {
   pMaterial  pm;

src/scene.c

@@ -1,14 +1,72 @@
 #include "medit.h"
 #include "extern.h"
 #include "sproto.h"
-
+#include "mesh.h"
 
 extern GLboolean  hasStereo;
 extern int       *pilmat,ipilmat,refmat,reftype,refitem;
 extern short      schw,schh;
 extern ubyte      quiet,fullscreen,tiling,stereoMode;
 
 
+void displayTetraMetrics(pScene sc) {
+    if (!sc->displayTetraMetrics)
+        return;
+
+    // Saving the projection matrix
+    glMatrixMode(GL_PROJECTION);
+    glPushMatrix();
+    glLoadIdentity();
+    // Configuring an orthographic matrix
+    gluOrtho2D(0, sc->par.xs, 0, sc->par.ys);
+
+    // Saving the matrix MODELVIEW
+    glMatrixMode(GL_MODELVIEW);
+    glPushMatrix();
+    glLoadIdentity();
+
+    // Choosing a color for the text (white here)
+    glColor3f(1.0, 1.0, 1.0);
+
+    // Position of the text : here (10, sc->par.ys - 20)
+    // sc->par.xs is the width of the window, sc->par.ys its height.
+    // So we place the text 10 pixels away from the left side & 20 pixels below the upper side.
+    glRasterPos2i(10, sc->par.ys - 20);
+
+    // Preparing the string for the tetrahedron quality
+    char info[256];
+    sprintf(info, "Tetra %d: Quality = %g", sc->currentTetraMetrics.id, sc->currentTetraMetrics.quality);
+    // Positioning the text
+    for (int i = 0; info[i] != '\0'; i++) {
+        glutBitmapCharacter(GLUT_BITMAP_HELVETICA_18, info[i]);
+    }
+
+    // Positioning a 2nd line to display the edge lengths.
+    glRasterPos2i(10, sc->par.ys - 40);
+    sprintf(info, "Edge lengths: %g %g %g %g %g %g",
+            sc->currentTetraMetrics.edgeLengths[0],
+            sc->currentTetraMetrics.edgeLengths[1],
+            sc->currentTetraMetrics.edgeLengths[2],
+            sc->currentTetraMetrics.edgeLengths[3],
+            sc->currentTetraMetrics.edgeLengths[4],
+            sc->currentTetraMetrics.edgeLengths[5]);
+    // Displaying the edge lengths
+    for (int i = 0; info[i] != '\0'; i++) {
+        glutBitmapCharacter(GLUT_BITMAP_HELVETICA_18, info[i]);
+    }
+
+    // Restoring the matrix MODELVIEW en pop
+    glPopMatrix();
+
+    // Restoring the projection matrix and pulling it out of 2D mode
+    glMatrixMode(GL_PROJECTION);
+    glPopMatrix();
+
+    // Coming back to MODELVIEW default mode
+    glMatrixMode(GL_MODELVIEW);
+}
+
+
 /* return current active scene */
 int currentScene() {
   int  k,idw;
@@ -651,6 +709,8 @@ void redrawScene() {
     glTranslatef(sc->cx,sc->cy,sc->cz);
     drawModel(sc);
     if ( sc->type & S_DECO )  redrawStatusBar(sc);
+
+    displayTetraMetrics(sc);
   }
 
   else {
@@ -686,6 +746,8 @@ void redrawScene() {
     drawModel(sc);
     if ( sc->type & S_DECO )  redrawStatusBar(sc);
 
+    displayTetraMetrics(sc);
+
     /* right view */
     glDrawBuffer(GL_BACK_RIGHT);
     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
@@ -707,8 +769,11 @@ void redrawScene() {
     glTranslatef(sc->cx,sc->cy,sc->cz);
     drawModel(sc);
     if ( sc->type & S_DECO )  redrawStatusBar(sc);
+
+    displayTetraMetrics(sc);
   }
 
+
   /* refresh screen */
   if ( saveimg && animate )
     glutSwapBuffers();

src/sproto.h

@@ -143,6 +143,8 @@ void   meshCoord(pMesh ,int );
 void   meshBox(pMesh mesh,int bb);
 void   meshRef(pScene sc,pMesh mesh);
 int    meshUpdate(pScene sc,pMesh mesh);
+TetraMetrics compute_mesh_quality_edges_length(pMesh mesh, int tetraIndex);
+
 
 /* mlists.c */
 GLuint  listTriaMap(pScene ,pMesh );
@@ -165,6 +167,7 @@ void   redrawOverlay(int stretchX,int stretchY);
 void   motionCamera(int x,int y);
 void   mouseCamera(int button,int state,int x,int y);
 void   animateCamera();
+extern int reftype;
 
 /* normal.c */
 GLuint drawNormals(pMesh mesh,pScene sc);
@@ -202,6 +205,7 @@ GLuint pickingList(pScene ,int ,int );
 GLuint pickingPoint(pScene sc,int x,int y);
 GLuint pickItem(pMesh ,pScene ,int );
 GLuint pickingScene(pScene sc,int x,int y,int ident);
+int pickingTetrahedron(pScene sc, int x, int y);
 
 /* prierr.c */
 void   prierr(int typerr,int indice);
@@ -226,6 +230,7 @@ void   deleteScene(pScene sc);
 void   initGrafix(pScene sc,pMesh mesh);
 int    createScene(pScene sc,int idmesh);
 void   streamIdle();
+void displayTetraMetrics(pScene sc);
 
 /* scissor.c */
 void   scissorScene();