Game Physics: Implementing Support Function for Polyhedrons using Half Edges

This post is part of my Game Physics Series.

The half edge data structure is a nice data structure to represent polyhedron meshes. I will first introduce to you the half edge data structure; then, I will show you an example implementation of the half edge data structure; lastly, I will explain how we can implement an efficient support function for polyhedrons using a hill-climbing technique with the half edge data structure.

Please note that this entire post assumes convex polyhedrons.

The Half Edge Data Structure

The half edge data structure is used to represent meshes using 3 lists: a list of vertices, a list of half edges (directed edges, distinguished from regular bi-directional edges), and a list of faces.

A vertex is a point in 3D. Each half edge points from one vertex to another; two half edges consisted of the same pair of vertices but pointing in opposite directions are twins. A face is a triangle. Vertices, edges, and faces are all features of a mesh.

Note that a face does not necessarily have to be a triangle in a half edge data structure; it can be a polygon of any number of vertices. However, for simplicity’s sake, the example implementation shown in this post uses triangle faces.

The implementation I will show below is index-based instead of pointer-based in terms of reference-keeping, because this implementation supports dynamic addition and removal of vertices and faces, managing its internal memory usage. A free list variable for a feature stores the index of a free feature, and each free feature has a free link property that refers to the next free feature; basically, a free list is implemented as a singly linked list.

We define a face by a loop of 3 half edges, where the face is to the left of each half edge (CCW winding).

Each half edge contains the index of the face to its left, the index of the vertex it is pointing at, the index of the next half edge of the same face, the index of the previous half edge of the same face, and the index of its twin. Each vertex contains a position and the index of an outgoing half edge; it doesn’t matter which half edge it is, for we can recover the entire set of outgoing edges. Each face contains the index of an edge that the face is formed by; same as half edges, it doesn’t matter which edge it is, because we can always recover the entire loop of half edges forming the face.

The half edge data structure essentially contains adjacency information of features. We can perform various adjacency queries on such structure, such as querying for all adjacent vertices of a vertex, or querying for the two adjacent faces of an edge. Hence, I use the name AdjacencyInfo for the data structure.

Below is the interface for the AdjacencyInfo structure.

struct AdjacencyInfo
{
  // half edge
  struct Edge
  {
    int vert; // pointed to by half edge
    int face; // to the left of half edge
    int next;
    int prev;
    int twin;
    int freeLink;
    bool active;
    Edge(void)
      : vert(-1)
      , face(-1)
      , next(-1)
      , prev(-1)
      , twin(-1)
      , freeLink(-1)
      , active(false)
    { }
  };

  struct Vert
  {
    Vec3 position;
    int edge; // a half edge pointing away
    int freeLink;
    bool active;
    Vert(void)
      : edge(-1)
      , freeLink(-1)
      , active(false)
    { }
  };

  struct Face
  {
    int edge; // an incident half edge
    int freeLink;
    bool active;
    Face(void) 
      : edge(-1)
      , freeLink(-1)
      , active(false)
    { }
  };

  // type defs
  typedef std::vector<Vert> VertList;
  typedef std::vector<Edge> EdgeList;
  typedef std::vector<Face> FaceList;

  //constructure
  AdjacencyInfo(void) { Clear(); }

  // operations
  int AddVert(const Vec3 &position);
  void RemoveVert(int v);
  int AddFace(int v0, int v1, int v2);
  void RemoveFace(int f);
  void Clear(void);

  // utils
  int FindVertEdge(int v) const;

  // container for memory management
  VertList m_verts;
  EdgeList m_edges;
  FaceList m_faces;

  // free lists
  int m_freeVerts;
  int m_freeEdges;
  int m_freeFaces;

  // counters
  uint m_numVerts;
  uint m_numEdges;
  uint m_numFaces;

  // edge records
  typedef std::pair<int, int> VertPair;
  typedef std::map<VertPair, int> EdgeMap;
  EdgeMap m_edgeMap;
};

First, let’s look at the function templates that allocate and free features. A feature’s active property is set in these two function templates, indicating whether this feature is allocated from the container or not.

template <typename T, typename Container>
int Allocate(int &freeList, 
             Container &container)
{
  int index = -1;
  if (freeList < 0)
  {
    // free list empty, just push
    index = static_cast<int>(container.size());
    container.resize(index + 1);
  }
  else
  {
    index = freeList;

    // grab free feature    
    T &feature = container[freeList]; 

    // free list
    freeList = feature.freeLink;

    // remove feature from free list
    feature.freeLink = -1;
  }
  T &feature = container[index];
  feature.active = true;
  return index;
}

template <typename T, typename Container>
static void Free(int index, 
                 int &freeList, 
                 Container &container)
{
  T &feature = container[index];
  feature.freeLink = freeList;
  feature.active = false;
  freeList = index;
}

Now, let’s look at the methods that adds a new vertex and removes a vertex from the mesh. Notice that the removal method destroys faces that are connected to the vertex.

int AdjacencyInfo::AddVert(const Vec3 &position)
{
  const int v = Allocate<Vert>(m_freeVerts, m_verts);
  auto &vert = m_verts[v];

  // write position
  vert.position = position;

  // initialize edge index
  vert.edge = -1;

  ++m_numVerts;

  return v;
}

void AdjacencyInfo::RemoveVert(int v)
{
  auto &vert = m_verts[v];

  // remove faces this vert connects to
  const int startEdge = vert.edge;
  if (startEdge >= 0)
  {
    int e = startEdge;
    do 
    {
      auto &edge = m_edges[e];
      RemoveFace(edge.face);
      e = m_edges[edge.prev].twin;
    } while (e != startEdge);
  }

  // dispose vert
  if (vert.edge >= 0)
    vert.edge = -1;

  // update free list
  Free<Vert>(v, m_freeVerts, m_verts);

  --m_numVerts;
}

Next, let’s look at the methods that adds a new face and removes a face from the mesh. The face creation requires 3 indices of vertices forming the face. This process is a little bit more complicated, so I will explain the algorithm with the comments.

int AdjacencyInfo::AddFace(int v0, int v1, int v2)
{
  const int faceVerts[3] = { v0, v1, v2 };
  
  // allocate face
  const int f = Allocate<Face>(m_freeFaces, m_faces);
  auto &face = m_faces[f];

  // create edges
  {
    // iterate face edges
    int faceEdgeIndices[3] = {-1, -1, -1};
    for (uint i = 2, j = 0; j < 3; i = j++)
    {
      const uint v0 = faceVerts[i];
      const uint v1 = faceVerts[j];
      auto &vert0 = m_verts[v0];
      auto &vert1 = m_verts[v1];

      // check existence of half edge pair
      const auto edgeIter = 
        m_edgeMap.find(VertPair(v0, v1));
      const bool edgePairExists = 
        edgeIter != m_edgeMap.end();
      int e01 = -1;
      int e10 = -1;
      if (edgePairExists)
      {
        e01 = edgeIter->second;
        e10 = m_edges[e01].twin;
      }
      else
      {
        // allocate & init half edge pair
        e01 = Allocate<Edge>(m_freeEdges, m_edges);
        e10 = Allocate<Edge>(m_freeEdges, m_edges);

        // link twins
        m_edges[e01].twin = e10;
        m_edges[e10].twin = e01;

        // record edge existence
        m_edgeMap[VertPair(v0, v1)] = e01;
        m_edgeMap[VertPair(v1, v0)] = e10;

        m_numEdges += 2;
        m_numBoundaryEdges += 2;
      } // end of edge allocation

      auto &edge01 = m_edges[e01];
      auto &edge10 = m_edges[e10];

      // link vert to edges
      if (edge01.vert < 0)
        edge01.vert = v1;
      if (edge10.vert < 0)
        edge10.vert = v0;

      // link face to edge
      if (edge01.face < 0)
      {
        edge01.face = f;
        --m_numBoundaryEdges;
      }
      
      // link edge to vert
      if (vert0.edge < 0)
        vert0.edge = e01;

      // link edge to face
      if (face.edge < 0)
        face.edge = e01;

      // record face edges
      faceEdgeIndices[i] = e01;
    }

    // link face edges
    for (unsigned i = 2, j = 0; j < 3; i = j++)
    {
      const int eI = faceEdgeIndices[i];
      const int eJ = faceEdgeIndices[j];
      m_edges[eI].next = eJ;
      m_edges[eJ].prev = eI;
    }

  } // end of edge creation

  ++m_numFaces;
  return f;
}

int AdjacencyInfo::FindVertEdge(int v) const
{
  // POSSIBLE OPTIMIZATION: 
  //   this is currently a linear operation
  for (auto &pair : m_edgeMap)
  {
    if (v == pair.first.first)
      return pair.second;
  }
  return -1;
}

void AdjacencyInfo::RemoveFace(int f)
{
  // remove adjacent edges
  const int startEdge = m_faces[f].edge;

  int faceVertIndices[3] = {-1, -1, -1};
  int faceEdgeIndices[3] = {-1, -1, -1};
  int e = startEdge;
  int i = 0;
  do
  {
    auto &edge = m_edges[e];
    const int t = edge.twin;
    const int n = edge.next;
    const int p = edge.prev;
    auto &twin = m_edges[t];
    auto &next = m_edges[n];
    auto &prev = m_edges[p];

    faceVertIndices[i] = edge.vert;
    faceEdgeIndices[i] = e;
    ++i;

    // free both edges if twin face does not exist
    if (twin.face < 0)
    {
      m_edgeMap.erase(VertPair(edge.vert, prev.vert));
      m_edgeMap.erase(VertPair(prev.vert, edge.vert));

      edge.twin = -1;
      twin.twin = -1;

      Free<Edge>(e, m_freeEdges, m_edges);
      Free<Edge>(t, m_freeEdges, m_edges);

      m_numEdges -= 2;
      m_numBoundaryEdges -= 2;
    }

    ++m_numBoundaryEdges;

    e = n;
  } while (e != startEdge);

  // unlink everything from edges
  for (const int e : faceEdgeIndices)
  {
    // okay to access data members after fake &quot;free&quot;
    auto &edge = m_edges[e];
    edge.next = -1;
    edge.prev = -1;
    edge.vert = -1;
    edge.face = -1;
  }

  // update vert edge
  for (const int v : faceVertIndices)
  {
    auto &vert = m_verts[v];
    vert.edge = FindVertEdge(v);
  }

  // unlink edge from face
  auto &face = m_faces[f];
  face.edge = -1;

  // finally, free face
  Free<Face>(f, m_freeFaces, m_faces);

  --m_numFaces;
}

Lastly, the trivial clear method that removes all features:

void AdjacencyInfo::Clear(void)
{
  m_verts.clear();
  m_edges.clear();
  m_faces.clear();

  m_numVerts = 0;
  m_numEdges = 0;
  m_numFaces = 0;

  m_freeVerts = -1;
  m_freeEdges = -1;
  m_freeFaces = -1;

  m_edgeMap.clear();
}

Adjacency Query

As mentioned before, we are ultimately going to use the half edge data structure to implement a support function for polyhedrons using a hill-climbing approach. Hill-climbing requires the ability to query adjacent vertices of a given vertex. Below is the general algorithm for performing an operation on all adjacent vertices of a vertex.

// loop through all edges leaving myVert
int e = myVert->edge;
int startEdge = e;
const float oldBestDot = bestDot;
do
{
  auto &edge = edges[e];

  // get adjacent vert pointed by edge
  const int vertIndex = edge.vert;
  auto vert = &verts[vertIndex];
  
  // PERFORM OPERATION ON vert HERE
  
  // get next out-going edge
  e = edges[edge.twin].next;
} while (startEdge != e);

Hill-Climbing Support Function

If you’re implementing a polyhedron collider with traditional triangle lists, then you’ll probably end up using some sort of brute-force approach to implement its support function.

const Vec3 Support(const Mesh &mesh, const Vec3 &dir)
{
  float bestDot = -FLT_MAX;
  Vec3 support;
  for (const auto &vert : mesh.Verts())
  {
    const float tempDot = vert.Dot(dir);
    if (tempDot > bestDot)
    {
      bestDot = tempDot;
      support = vert;
    }
  }
  return support;
}

If the polyhedron collider is stored as a half edge structure, then we can use a hill-climbing technique to implement the support function.

The concept is simple: We start with a random vertex (the first vertex in the vertex list will do) and compute its dot product with the support direction (support value). We then use the adjacency information to find all the vertices around the vertex, compute the support values of these vertices, move to the vertex that has the largest support value, and repeat this process until no adjacent vertex of the current vertex can give a better support value. The last vertex is our support point.

const Vec3 Support(const Mesh &mesh, const Vec3 &dir)
{
  // grab adjacency information
  auto &adj = mesh.GetAdjacencyInfo();
  auto &verts = adj.m_verts;
  auto &edges = adj.m_edges;
  
  // start from first point
  int bestVertIndex = 0;
  auto bestVert = &verts[bestVertIndex];
  float bestDot = dir.Dot(bestVert->position);

  bool end = false;
  while (!end)
  {
    // loop through all edges leaving best vert
    int e = bestVert->edge;
    int startEdge = e;
    const float oldBestDot = bestDot;
    do
    {
      auto &edge = edges[e];

      // get adjacent vert pointed by edge
      const int vertIndex = edge.vert;
      auto vert = &verts[vertIndex];

      const float dot = dir.Dot(vert->position);
      if (dot >= bestDot + EPSILON)
      {
        // record new best vert
        bestVertIndex = bestVertIndex;
        bestVert = vert;
        bestDot = dot;
      }

      // get next out-going edge
      e = edges[edge.twin].next;
    } while (startEdge != e);

    // end if no new best vert is found
    end = (bestDot == oldBestDot);
  }

  return bestVert->position;
}

End of Half Edge Data Structure for Polyhedrons

Now that you understand the half edge data structure, you should be able to implement a support function for polyhedrons that performs better than the naive brute-force approach (with reasonably many vertices per polyhedron, that is).

About Allen Chou

Physics / Graphics / Procedural Animation / Visuals
This entry was posted in Gamedev, Physics. Bookmark the permalink.

Leave a Reply