Search code examples
c#wpfcomputational-geometryconvexconcave-hull

Concave hull algorithm from pseudocode to C#


I am trying to convert the algorithm as described here (page 12) from pseudocode into working C# code. The algorithm describes how a convex hull is 'transformed' into a concave hull by breaking up edges that are considered too long into smaller edges. I understand the general idea that the authors propose, but have trouble converting this into working code. Please see below the code that I have got so far, including comments (//) at the start of each pseudocode line. The problem I am having is not so much with a specific line - though I am certain that the current way of calculating 'localMaximumDistance' is not correct. If anyone has any pointers on how to approach this I would really like to hear those. (In pseudocode, this is the line that says 'calculate local maximum distance d for edges;')

Thank you in advance for your time and feedback! :)

List<Line> concaveLineList = new List<Line>();
List<Line> sortedList = lineList.OrderByDescending(CalculateLength).ToList();
const int concaveTreshhold = 40;

PointCollection concavePointCollection = new PointCollection();
while (sortedList.Count > 0) {
    Console.WriteLine(concaveLineList.Count.ToString());
    //select the longest edge e from list A
    Line longestLine = sortedList[0];
    //remove edge e from list A
    sortedList.RemoveAt(0);
    //calculate local maximum distance d for edges - ???
    //double localMaximumDistance = CalculateLength(longestLine);
    List<Point<double>> nearbyPoints = new List<Point<double>>();

    foreach (BallUc ballUc in ballUcList) {
        if (Math.Abs(ballUc .CurrentPosition.X - longestLine.X1) > concaveTreshhold &&
        Math.Abs(ballUc .CurrentPosition.Y - longestLine.Y1) > concaveTreshhold) { 
            nearbyPoints.Add(new Point<double>(ballUc.CurrentPosition.X, ballUc.CurrentPosition.Y));
            }
        }

        double lineLenght = CalculateLength(longestLine);
        double localMaximumDistance = lineLenght / nearbyPoints.Count + concaveTreshhold * 4; //this value is based on nothing currently..

        if (lineLenght > localMaximumDistance) {
            //find the point p with the smallest maximum angle a
            Point smallestAnglePoint = new Point();
            double? smallestAngle = null;
            foreach (Point p in pointCollection) {
                if ((p.X == longestLine.X1 && p.Y == longestLine.Y1) ||
                (p.X == longestLine.X2 && p.Y == longestLine.Y2)) {
                    //these are the points already in the line.
                }
                else {
                    Line tempLine1 = new Line {X1 = p.X, X2 = longestLine.X1, Y1 = p.Y, Y2 = longestLine.Y1};
                    Line tempLine2 = new Line {X1 = p.X, X2 = longestLine.X2, Y1 = p.Y, Y2 = longestLine.Y2};

                    //calculate angle between the longest edge and the new edges
                    double angleInRadians1 = Math.Atan2(p.Y, p.X) - Math.Atan2(tempLine1.Y2, tempLine1.X2);
                    double angleInRadians2 = Math.Atan2(p.Y, p.X) - Math.Atan2(tempLine2.Y2, tempLine2.X2);
                    //select the largest angle of the two angles
                    double largestLocalAngle = Math.Max(angleInRadians1, angleInRadians2);

                    //in case of first calculation, smallestAngle is still null - in this case it should be assigned the value
                    //(this is probably not very elegant code)
                    if (smallestAngle == null) {
                        smallestAngle = largestLocalAngle;
                        smallestAnglePoint = p;
                    }
                    //we have to find the smallest angle.
                    else if (largestLocalAngle < smallestAngle) {
                        smallestAngle = largestLocalAngle;
                        smallestAnglePoint = p;
                    }
                    //double angleinDegrees = angleInRadians * 180 / Math.PI;
                }
            }
            //TODO if angle a is small enough and point p is not on the boundary

            //create edges e2 and e3 between point p and endpoints of edge e
            Line e2 = new Line {
                X1 = smallestAnglePoint.X,
                Y1 = smallestAnglePoint.Y,
                X2 = longestLine.X1,
                Y2 = longestLine.Y1
            };
            sortedList.Add(e2);
            Line e3 = new Line {
                X1 = smallestAnglePoint.X,
                Y1 = smallestAnglePoint.Y,
                X2 = longestLine.X2,
                Y2 = longestLine.Y2
            };
            sortedList.Add(e3);

            //if edge e2 and e3 don't intersect any other edge
            foreach (Line line in sortedList) {
                Point lineInitialPoint = new Point(line.X1, line.Y1);
                Point lineTerminalPoint = new Point(line.X2, line.Y2);

                Point line2InitialPoint = new Point(e2.X1, e2.Y1);
                Point line2TerminalPoint = new Point(e2.X2, e2.Y2);

                Point line3InitialPoint = new Point(e2.X1, e2.Y1);
                Point line3TerminalPoint = new Point(e2.X2, e2.Y2);

                Point intersectionPoint = GetIntersection(line2InitialPoint, line2TerminalPoint, lineInitialPoint, lineTerminalPoint);
                Point intersectionPoint2 = GetIntersection(line3InitialPoint, line3TerminalPoint, lineInitialPoint, lineTerminalPoint);

                 if ((Double.IsNaN(intersectionPoint.X) && Double.IsNaN(intersectionPoint.Y)) &&
                    (Double.IsNaN(intersectionPoint2.X) && Double.IsNaN(intersectionPoint2.Y))) {
                     //no intersection found, keep rolling..
                     //Console.WriteLine("no intersection found");
                 }
                 else {
                    //intersection found, lines no longer valid
                    Console.WriteLine("intersection found");
                    break;
                 }
                 concaveLineList.Add(e2);
                 concaveLineList.Add(e3);
            }
        }
        //if edge e2 and e3 was not added to list A
        else {
            //add edge e to list B
            concaveLineList.Add(longestLine);
            concavePointCollection.Add(new Point(longestLine.X1, longestLine.Y1));
            concavePointCollection.Add(new Point(longestLine.X2, longestLine.Y2));
        }
    }

Solution

  • I have implemented this algorithm via JavaScript, may be it will help you. There you can see function that transforms convex hull to concave hull. My implementation of the algorithm is not exactly the same as in the Paper, but it based on it.

    In my implementation I skipped localMaximumDistance calculation :-) and decided that first of all I need to see use case in which my algorithm without localMaximumDistance will work wrong.