Volko Hull - An Alternative to Convex Hull

When you have a point cloud and want to obtain the shape, you can use the convex hull or the alpha shape algorithm. In the computational geometry library I am using in my projects, CGAL, these algorithms are built-in, but they have a disadvantage: the points are not output in the correct order. I first tried to solve this problem by implementing a convex hull algorithm myself that outputs the points in the correct order, but this did not always help since it ignored concave corners. Therefore I wrote my own algorithm to replace it.

My algorithm divides the rectangular area of the given cluster into four triangles and computes the points belonging to the outer border of each triangle. Thus it also includes concave corners.

Imagine this to be the rectangle:

--------
|\    /|
| \ A/ |
| C\/  |
|  /\D |
| /B \ |
|/    \|
--------

If it is a square, we can assume that its sides have a length of 2a. The coordinates of the mid point are always (0, 0).

Then, the points have to fulfil the following conditions to belong to one of the four triangles:

A: |x| + y <= a
B: |x| - y <= a
C: -x + |y| <= a
D: +x + |y| <= a

Generalizing this for a rectangle with horizontal length 2a and vertical length 2b, this becomes:

A: |x| + y a/b <= a
B: |x| - y a/b <= a
C: -x b/a + |y| <= b
D: +x b/a + |y| <= b

Now you should be able to understand the main method:

// computes the Volko hull of a cluster
std::vector<Point_3> compute_volko_hull(std::vector<Point_3> cluster)
{
  std::vector<Point_2> cluster_full = convert_cluster(cluster);
  auto bbox = CGAL::bounding_box(cluster_full.begin(), cluster_full.end());
  double dx = (bbox.xmax() - bbox.xmin()) / 2;
  double dy = (bbox.ymax() - bbox.ymin()) / 2;
  Point_2 midPoint = Point_2(bbox.xmin() + dx, bbox.ymin() + dy);

  std::vector<Point_2> cluster_A;
  std::vector<Point_2> cluster_B;
  std::vector<Point_2> cluster_C;
  std::vector<Point_2> cluster_D;

  for (const auto& pt : cluster_full)
  {
    double lx = pt.x() - midPoint.x();
    double ly = pt.y() - midPoint.y();

    if (abs(lx) + ly * dx / dy <= dx)
      cluster_A.push_back(pt);
    else if (abs(lx) - ly * dx / dy <= dx)
      cluster_B.push_back(pt);
    else if (-lx * dy / dx + abs(ly) <= dy)
      cluster_C.push_back(pt);
    else if (+lx * dy / dx + abs(ly) <= dy)
      cluster_D.push_back(pt);
  }

  std::sort(cluster_A.begin(), cluster_A.end());
  cluster_A = removedupes(cluster_A, false);
  std::sort(cluster_B.rbegin(), cluster_B.rend()); // reverse
  cluster_B = removedupes(cluster_B, true);
  cluster_C = swapxy(cluster_C);
  std::sort(cluster_C.rbegin(), cluster_C.rend()); // reverse
  cluster_C = removedupes(cluster_C, false);
  cluster_C = swapxy(cluster_C);
  cluster_D = swapxy(cluster_D);
  std::sort(cluster_D.begin(), cluster_D.end());
  cluster_D = removedupes(cluster_D, true);
  cluster_D = swapxy(cluster_D);

  std::vector<Point_2> cluster_result;
  for (const auto& c : cluster_A)
    cluster_result.push_back(c);
  for (const auto& c : cluster_D)
    cluster_result.push_back(c);
  for (const auto& c : cluster_B)
    cluster_result.push_back(c);
  for (const auto& c : cluster_C)
    cluster_result.push_back(c);

  return convert_cluster_back(cluster_result);
}

We need two auxiliary methods, one for swapping the x and y coordinates (because std::sort prioritizes x and we sometimes need to prioritize y) and one for removing the inner points. These can be implemented as follows:

// auxiliary method for compute_volko_hull
std::vector<Point_2> swapxy(std::vector<Point_2> cluster)
{
  std::vector<Point_2> newCluster;
  for (const auto &pt : cluster) newCluster.push_back(Point_2(pt.y(), pt.x()));
  return newCluster;
}

// auxiliary method for compute_volko_hull
std::vector<Point_2> removedupes(std::vector<Point_2> cluster, bool chooseMax)
{
  if (cluster.size() < 2) return cluster;

  std::vector<Point_2> newCluster;
  for (int i = 0; i < cluster.size(); )
  {
    double yMin = cluster[i].y();
    double yMax = cluster[i].y();
    int j = i + 1;
    for (; j < cluster.size(); j++)
    {
      if (abs(cluster[j].x() - cluster[i].x()) < .00001)
      {
        if (cluster[j].y() < cluster[i].y())
        {
          if (cluster[j].y() < yMin)
            yMin = cluster[j].y();
        }
        else
        {
          if (cluster[j].y() > yMax)
            yMax = cluster[j].y();
        }
      }
      else
        break;
    }

    if (chooseMax)
      newCluster.push_back(Point_2(cluster[i].x(), yMax));
    else
      newCluster.push_back(Point_2(cluster[i].x(), yMin));

    i = j;
  }
  return newCluster;
}

To compile this code directly, you will need a reasonably recent version of CGAL.