Implement convex hull algorithm.

README.md

@@ -36,7 +36,7 @@ Two caps and many strips composites a cylinder.
 Almost all 3D editor have a infinite grid ground, I just made a finite one, in the future, I should expand the grid outside of the screen to make it infinite.  
 Now, for just beginning, I think it's a not bad start.
 <img src="screenshot/dust3d_sphere_cylinder.png">
-- [ ] Implement B-Mesh algorithm (Dec 18, 2016 ~ Dec 25, 2016)   
+- [ ] Implement B-Mesh algorithm (Dec 18, 2016 ~ Dec 26, 2016)   
 *Drawing Skeletal Shape Balls*  
 Draw shape ball is easy, no need to rotate, I just need scale it along the ball's radius.
 Draw the cylinder which connects two shape balls is more difficult, I need do some math to rotate it. [Here](http://www.thjsmith.com/40/cylinder-between-two-points-opengl-c) described it.  
@@ -54,8 +54,8 @@ I created the test nodes's geometry information from Blender. Here is the render
 <img src="screenshot/dust3d_generate_quad.png" width="124" height="128">  
 When I am implementing the B-Mesh algorithm, I am also think in the future, how to create a library of bunch of initial base models. There is a paper [the Skeleton of a Closed 3D Shape](http://www1.idc.ac.il/icgf/GraphicsSeminar2006/DCGskeleton06.pdf) described how to generate skeleton from mesh, this is the reverse progress of what I am doing, I think it can resolve the problem of insufficient initial base models, I can generate from tons of existed models.  
 *Convex Hull*  
-After finish the rotation at the two connected bones, I need implement [3D Convex Hull](https://en.wikipedia.org/wiki/Convex_hull) algorithm at joint ball.
+After finish the rotation at the two connected bones, I need implement 3D Convex Hull algorithm at the joint ball, there are so many methods to get the convex hull, I found the [Gift wrapping](http://dccg.upc.edu/people/vera/wp-content/uploads/2014/11/GA2014-ConvexHulls3D-Roger-Hernando.pdf) is the most strait-forward one, though is not the most efficient one.  
-
+<img src="screenshot/dust3d_convex_hull.png" width="124" height="128">
 - [ ] Export Wavefront .obj  
 - [ ] Render B-Mesh result  
 - [ ] Design UI for monster parts configuration  

build/dust3d.pro

@@ -14,7 +14,8 @@ SOURCES += main.cpp \
   draw.cpp \
   array.c \
   bmesh.c \
-  matrix.c
+  matrix.c \
+  convexhull.c
 
 HEADERS += mainwindow.h \
 	render.h \
@@ -22,4 +23,5 @@ HEADERS += mainwindow.h \
   draw.h \
   array.h \
   bmesh.h \
-  matrix.h
+  matrix.h \
+  convexhull.h

src/bmesh.c

@@ -6,6 +6,7 @@
 #include "bmesh.h"
 #include "array.h"
 #include "matrix.h"
+#include "convexhull.h"
 #include "draw.h"
 
 #define BMESH_STEP_DISTANCE 0.4
@@ -245,6 +246,7 @@ static int bmeshGenerateInbetweenBallsBetween(bmesh *bm,
   generateYZfromBoneDirection(&boneDirection,
     &localYaxis, &localZaxis);
   
+  /*
   glColor3f(0.0, 0.0, 0.0);
   drawDebugPrintf("<%f,%f,%f> <%f,%f,%f> <%f,%f,%f>",
     localYaxis.x,
@@ -256,6 +258,7 @@ static int bmeshGenerateInbetweenBallsBetween(bmesh *bm,
     boneDirection.x,
     boneDirection.y,
     boneDirection.z);
+  */
   
   distance = vec3Length(&boneDirection);
   if (distance > BMESH_STEP_DISTANCE) {
@@ -485,7 +488,7 @@ static bmeshBall *bmeshFindBallForConvexHull(bmesh *bm, bmeshBall *root,
   bmeshBallIterator iterator;
   bmeshBall *child;
   float distance = vec3Distance(&root->position, &ball->position);
-  if (distance >= root->radius) {
+  if (distance > root->radius) {
     return ball;
   }
   child = bmeshGetBallFirstChild(bm, ball, &iterator);
@@ -502,11 +505,51 @@ static int bmeshStichFrom(bmesh *bm, bmeshBall *ball) {
   bmeshBall *child;
   bmeshBall *ballForConvexHull;
   if (BMESH_BALL_TYPE_ROOT == ball->type) {
+    convexHull *hull = convexHullCreate();
+    if (!hull) {
+      fprintf(stderr, "%s:convexHullCreate failed.\n", __FUNCTION__);
+      return -1;
+    }
     for (child = bmeshGetBallFirstChild(bm, ball, &iterator);
         child;
         child = bmeshGetBallNextChild(bm, ball, &iterator)) {
+      vec3 z, y;
+      quad q;
+      int vertexIndices[4];
+      
       ballForConvexHull = bmeshFindBallForConvexHull(bm, ball, child);
+      
+      vec3Scale(&ballForConvexHull->localYaxis, ballForConvexHull->radius, &y);
+      vec3Scale(&ballForConvexHull->localZaxis, ballForConvexHull->radius, &z);
+      vec3Sub(&ballForConvexHull->position, &y, &q.pt[0]);
+      vec3Add(&q.pt[0], &z, &q.pt[0]);
+      vec3Sub(&ballForConvexHull->position, &y, &q.pt[1]);
+      vec3Sub(&q.pt[1], &z, &q.pt[1]);
+      vec3Add(&ballForConvexHull->position, &y, &q.pt[2]);
+      vec3Sub(&q.pt[2], &z, &q.pt[2]);
+      vec3Add(&ballForConvexHull->position, &y, &q.pt[3]);
+      vec3Add(&q.pt[3], &z, &q.pt[3]);
+      
+      vertexIndices[0] = convexHullAddVertex(hull, &q.pt[0]);
+      vertexIndices[1] = convexHullAddVertex(hull, &q.pt[1]);
+      vertexIndices[2] = convexHullAddVertex(hull, &q.pt[2]);
+      vertexIndices[3] = convexHullAddVertex(hull, &q.pt[3]);
     }
+    convexHullGenerate(hull);
+    
+    glPushMatrix();
+    {
+      int triIndex;
+      glColor3f(1.0, 1.0, 1.0);
+      for (triIndex = 0; triIndex < convexHullGetTriangleNum(hull);
+          ++triIndex) {
+        triangle *tri = (triangle *)convexHullGetTriangle(hull, triIndex);
+        drawTriangle(tri);
+      }
+    }
+    glPopMatrix();
+    
+    convexHullDestroy(hull);
   }
   for (child = bmeshGetBallFirstChild(bm, ball, &iterator);
       child;

src/convexhull.c

@@ -0,0 +1,238 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include "convexhull.h"
+#include "array.h"
+#include "draw.h"
+
+//
+// Implement Gift wrapping method which describled in http://dccg.upc.edu/people/vera/wp-content/uploads/2014/11/GA2014-ConvexHulls3D-Roger-Hernando.pdf
+//
+//
+// Translate from Danielhst's lua version https://github.com/danielhst/3d-Hull-gift-wrap/blob/master/giftWrap.lua
+//
+
+struct convexHull {
+  array *vertexArray;
+  array *openEdgeArray;
+  array *triangleArray;
+  int nextEdgeIndex;
+  unsigned int *openEdgeExistMap;
+};
+
+typedef struct {
+  int firstVertex;
+  int secondVertex;
+} edge;
+
+convexHull *convexHullCreate(void) {
+  convexHull *hull = (convexHull *)calloc(1, sizeof(convexHull));
+  if (!hull) {
+    fprintf(stderr, "%s:Insufficient memory.\n", __FUNCTION__);
+    return 0;
+  }
+  hull->vertexArray = arrayCreate(sizeof(vec3));
+  if (!hull->vertexArray) {
+    fprintf(stderr, "%s:arrayCreate failed.\n", __FUNCTION__);
+    return 0;
+  }
+  hull->openEdgeArray = arrayCreate(sizeof(edge));
+  if (!hull->openEdgeArray) {
+    fprintf(stderr, "%s:arrayCreate failed.\n", __FUNCTION__);
+    return 0;
+  }
+  hull->triangleArray = arrayCreate(sizeof(triangle));
+  if (!hull->triangleArray) {
+    fprintf(stderr, "%s:arrayCreate failed.\n", __FUNCTION__);
+    return 0;
+  }
+  return hull;
+}
+
+void convexHullMarkEdgeAsExsits(convexHull *hull, int firstVertex,
+    int secondVertex) {
+  int mapIndex = firstVertex * arrayGetLength(hull->vertexArray) + secondVertex;
+  int unitIndex = mapIndex / (sizeof(unsigned int) * 8);
+  int unitOffset = mapIndex % (sizeof(unsigned int) * 8);
+  hull->openEdgeExistMap[unitIndex] |= (0x00000001 << unitOffset);
+}
+
+int convexHullEdgeExsits(convexHull *hull, int firstVertex,
+    int secondVertex) {
+  int mapIndex = firstVertex * arrayGetLength(hull->vertexArray) + secondVertex;
+  int unitIndex = mapIndex / (sizeof(unsigned int) * 8);
+  int unitOffset = mapIndex % (sizeof(unsigned int) * 8);
+  return hull->openEdgeExistMap[unitIndex] & (0x00000001 << unitOffset);
+}
+
+int convexHullAddVertex(convexHull *hull, vec3 *vertex) {
+  int newVertex = arrayGetLength(hull->vertexArray);
+  if (0 != arraySetLength(hull->vertexArray, newVertex + 1)) {
+    fprintf(stderr, "%s:arraySetLength failed.\n", __FUNCTION__);
+    return -1;
+  }
+  *((vec3 *)arrayGetItem(hull->vertexArray, newVertex)) = *vertex;
+  return newVertex;
+}
+
+int convexHullAddOpenEdgeNoCheck(convexHull *hull, int firstVertex,
+    int secondVertex) {
+  edge *e;
+  int newEdge = arrayGetLength(hull->openEdgeArray);
+  if (firstVertex < 0 || secondVertex < 0) {
+    fprintf(stderr, "%s:Invalid params(firstVertex:%d secondVertex:%d).\n",
+      __FUNCTION__, firstVertex, secondVertex);
+    return -1;
+  }
+  if (0 != arraySetLength(hull->openEdgeArray, newEdge + 1)) {
+    fprintf(stderr, "%s:arraySetLength failed.\n", __FUNCTION__);
+    return -1;
+  }
+  e = (edge *)arrayGetItem(hull->openEdgeArray, newEdge);
+  memset(e, 0, sizeof(edge));
+  e->firstVertex = firstVertex;
+  e->secondVertex = secondVertex;
+  return 0;
+}
+
+int convexHullAddEdge(convexHull *hull, int p1, int p2) {
+  convexHullMarkEdgeAsExsits(hull, p1, p2);
+  if (!convexHullEdgeExsits(hull, p2, p1)) {
+    return convexHullAddOpenEdgeNoCheck(hull, p2, p1);
+  }
+  return 0;
+}
+
+int convexHullAddTriangle(convexHull *hull, int firstVertex, int secondVertex,
+    int thirdVertex) {
+  triangle *tri;
+  int newTri = arrayGetLength(hull->triangleArray);
+  if (0 != arraySetLength(hull->triangleArray, newTri + 1)) {
+    fprintf(stderr, "%s:arraySetLength failed.\n", __FUNCTION__);
+    return -1;
+  }
+  tri = (triangle *)arrayGetItem(hull->triangleArray, newTri);
+  memset(tri, 0, sizeof(triangle));
+  tri->pt[0] = *((vec3 *)arrayGetItem(hull->vertexArray, firstVertex));
+  tri->pt[1] = *((vec3 *)arrayGetItem(hull->vertexArray, secondVertex));
+  tri->pt[2] = *((vec3 *)arrayGetItem(hull->vertexArray, thirdVertex));
+  return 0;
+}
+
+static void convexHullReleaseForGenerate(convexHull *hull) {
+  free(hull->openEdgeExistMap);
+  hull->openEdgeExistMap = 0;
+}
+
+void convexHullDestroy(convexHull *hull) {
+  arrayDestroy(hull->vertexArray);
+  arrayDestroy(hull->openEdgeArray);
+  arrayDestroy(hull->triangleArray);
+  convexHullReleaseForGenerate(hull);
+  free(hull);
+}
+
+static int convexHullPrepareForGenerate(convexHull *hull) {
+  free(hull->openEdgeExistMap);
+  hull->openEdgeExistMap = (unsigned int *)calloc(
+    arrayGetLength(hull->vertexArray) * arrayGetLength(hull->vertexArray) /
+      (sizeof(unsigned int) * 8) + 1,
+    sizeof(unsigned int));
+  if (!hull->openEdgeExistMap) {
+    fprintf(stderr, "%s:Insufficient memory.\n", __FUNCTION__);
+    return -1;
+  }
+  hull->nextEdgeIndex = 0;
+  return 0;
+}
+
+int convexHullGetLower(convexHull *hull) {
+  int index = 0;
+  int i;
+  for (i = 1; i < arrayGetLength(hull->vertexArray); ++i) {
+    vec3 *pI = (vec3 *)arrayGetItem(hull->vertexArray, i);
+    vec3 *pIndex = (vec3 *)arrayGetItem(hull->vertexArray, index);
+    if (pI->z < pIndex->z) {
+      index = i;
+    } else if (pI->z == pIndex->z) {
+      if (pI->y < pIndex->y) {
+        index = i;
+      } else if (pI->x < pIndex->x) {
+        index = i;
+      }
+    }
+  }
+  return index;
+}
+
+int convexHullGetNextVertex(convexHull *hull, int p1Index, int p2Index) {
+  vec3 pCorner = {1, 1, 0};
+  vec3 *p1 = (vec3 *)arrayGetItem(hull->vertexArray, p1Index);
+  vec3 *p2 = p2Index < 0 ? &pCorner : (vec3 *)arrayGetItem(hull->vertexArray,
+    p2Index);
+  vec3 edge;
+  int candidateIndex = -1;
+  int i;
+  vec3Sub(p2, p1, &edge);
+  vec3Normalize(&edge);
+  for (i = 0; i < arrayGetLength(hull->vertexArray); ++i) {
+    if (i != p1Index && i != p2Index) {
+      if (-1 == candidateIndex) {
+        candidateIndex = 0;
+      } else {
+        vec3 v, proj, candidate, canProj, cross;
+        vec3Sub((vec3 *)arrayGetItem(hull->vertexArray, i), p1, &v);
+        vec3ProjectOver(&v, &edge, &proj);
+        vec3Sub(&v, &proj, &v);
+        vec3Sub((vec3 *)arrayGetItem(hull->vertexArray,
+          candidateIndex), p1, &candidate);
+        vec3ProjectOver(&candidate, &edge, &canProj);
+        vec3Sub(&candidate, &canProj, &candidate);
+        vec3CrossProduct(&candidate, &v, &cross);
+        if (vec3DotProduct(&cross, &edge) > 0) {
+          candidateIndex = i;
+        }
+      }
+    }
+  }
+  return candidateIndex;
+}
+
+int convexHullGenerate(convexHull *hull) {
+  int index1, index2, index3;
+  convexHullReleaseForGenerate(hull);
+  if (0 != convexHullPrepareForGenerate(hull)) {
+    fprintf(stderr, "%s:convexHullPrepareForGenerate failed.\n", __FUNCTION__);
+    return -1;
+  }
+  index1 = convexHullGetLower(hull);
+  index2 = convexHullGetNextVertex(hull, index1, -1);
+  convexHullAddEdge(hull, index2, index1);
+  //glColor3f(0.0, 0.0, 0.0);
+  //drawDebugPrintf("edgeLength:%d", arrayGetLength(hull->openEdgeArray));
+  while (hull->nextEdgeIndex < arrayGetLength(hull->openEdgeArray)) {
+    edge *e = (edge *)arrayGetItem(hull->openEdgeArray, hull->nextEdgeIndex++);
+    index1 = e->firstVertex;
+    index2 = e->secondVertex;
+    if (convexHullEdgeExsits(hull, index1, index2)) {
+      continue;
+    }
+    convexHullMarkEdgeAsExsits(hull, index1, index2);
+    index3 = convexHullGetNextVertex(hull, index1, index2);
+    //drawDebugPrintf("%d,%d,%d", index1, index2, index3);
+    convexHullAddTriangle(hull, index1, index2, index3);
+    convexHullAddEdge(hull, index1, index2);
+    convexHullAddEdge(hull, index2, index3);
+    convexHullAddEdge(hull, index3, index1);
+  }
+  return 0;
+}
+
+int convexHullGetTriangleNum(convexHull *hull) {
+  return arrayGetLength(hull->triangleArray);
+}
+
+triangle *convexHullGetTriangle(convexHull *hull, int index) {
+  triangle *tri = (triangle *)arrayGetItem(hull->triangleArray, index);
+  return tri;
+}

src/convexhull.h

@@ -0,0 +1,14 @@
+#ifndef CONVEX_HULL_H
+#define CONVEX_HULL_H
+#include "vector3d.h"
+#include "draw.h"
+
+typedef struct convexHull convexHull;
+convexHull *convexHullCreate(void);
+int convexHullAddVertex(convexHull *hull, vec3 *vertex);
+void convexHullDestroy(convexHull *hull);
+int convexHullGenerate(convexHull *hull);
+int convexHullGetTriangleNum(convexHull *hull);
+triangle *convexHullGetTriangle(convexHull *hull, int index);
+
+#endif

src/vector3d.c

@@ -105,23 +105,27 @@ float vec3Angle(vec3 *a, vec3 *b) {
   float angle;
   vec3 p;
   float distance;
-  float dot;
   float acosParam;
-  float acosVal;
   vec3Sub(a, b, &p);
   distance = vec3Length(&p);
   if (0 == distance) {
     return 0;
   }
-  dot = vec3DotProduct(a, b);
+  acosParam = vec3DotProduct(a, b) / distance;
-  acosParam = dot / distance;
   if (acosParam < -1) {
     acosParam = -1;
   }
   if (acosParam > 1) {
     acosParam = 1;
   }
-  acosVal = acos(acosParam);
+  angle = 180 / M_PI * acos(acosParam);
-  angle = 180 / M_PI * acosVal;
   return angle;
 }
+
+void vec3ProjectOver(vec3 *a, vec3 *over, vec3 *result) {
+  float length = vec3DotProduct(a, over);
+  result->x = length * over->x;
+  result->y = length * over->y;
+  result->z = length * over->z;
+}
+

src/vector3d.h

@@ -24,6 +24,7 @@ float vec3Distance(vec3 *a, vec3 *b);
 void vec3Normal(vec3 *a, vec3 *b, vec3 *c, vec3 *normal);
 void vec3RotateAlong(vec3 *a, float angle, vec3 *axis, vec3 *result);
 float vec3Angle(vec3 *a, vec3 *b);
+void vec3ProjectOver(vec3 *a, vec3 *over, vec3 *result);
 
 #ifdef __cplusplus
 }