get nearest points that form a triangle

Question

I have some points (blue) in 2D.

I want to get these three points, that form a triangle in such a way, that the point D (red) lies inside this triangle. If there is no such triangle, an exception can be thrown.

So for the above picture I want to get the black points:

What I have done so far: I thought I could order the points by their distance to D, and than take the first three points out of that sorted list. But the problem is, that it could be, that these three nearest points form a triangle which does not include the point D. This is shown in this picture:

Additional to getting the wrong points, I could not determine weather or not D lies in the convex hull of the points found, and therefore if there is a triangle which includes the point D. This is the point, where I got stuck.

Answer 1

As correctly noted in the comments by TaW, the following algorithm in its basic form won't always find the best solution or a solution at all, because it starts greedily with the two closets points.

But this could be fixed with repeating the algorithm: if it won't find the triangle you can reiterate it ignoring the first closest point.

If you don't have many points you can repeat the algorithm always anyway for different starting points to be sure you'll find the optimal solution.

1) find the closest point to D, let's call it A

2) find the second closest point to D, let's call it B

3) find the equation of a line passing through D and A, let's call it L1 (the missing point must be on the other side of L1 than B)

4) find the equation of a line passing through D and B, let's call it L2 (the missing point must be on the other side of L2 than A)

5) filter the rest of points: leave only points that are on the other side of L1 than B, and that are on the other side of L2 than A

6) if there are no such points, throw exception (or reiterate with different starting point)

7) otherwise find the closest one, let's call it C

8) result is the triangle ABC

Additional notes:

Two points are on different sides of a line if this equation, given their coordinates gives different signs (X, Y, Z are line equation coefficients, usually A, B, C are used but I didn't want to confuse them with point labels above):

Equation of a line that passes through two points with coordinates (V1x, V1y) and (V2x, V2y) can be found with this formula:

Which gives you the following formulas for line equation coefficients:

Answer 2

The key is to look for valid solution while in the same time you focus on the optimal solution:

For each point: - store distance to target point - store position relative to target point.

I would use the following enum to store position:

enum RelativePosition
{
    ll,
    le,
    lg,
    eg,
    gg,
    ge,
    gl,
    el,
    ee
}

where first letter represents point's x coordinate relative to target's x coordinate and the second letter represents point's y coordinate relative to target's y coordinate.

l less, g greater, e equal

Order your points by distance (ascending) to target point Start with closest point and depending on the relative position get a list of candidates that would form a triangle around your target. Also, get closest to target from those candidates and continue with third point in the same manner.

I am on mobile now and it's hard to provide code but I'll be able to write it in one hour or two.

Edit:

Sorry for delay. Here is some code. You'll see that in ValidPositions method I hardcoded all the valid positions relative to first and second point's positions. I know that there is a mathematical relation between them and they can be generated but let's say I let that as an exercice. :) Even with this method there are cases when you can't be sure if the target point is in the area of the triangle (see UncertainSolution method). However the number of tests if TriangleContainsPoint is reduced.

Edit 2: Fixed bug in method TriangleContainsPoint

class Point2D
{
    public double X { get; set; }
    public double Y { get; set; }
}

enum RelPos2D
{
    ll = 1,
    le = 2,
    lg = 3,
    eg = 4,
    gg = 5,
    ge = 6,
    gl = 7,
    el = 8,
    ee = 0
}

static class Tools2D
{
    public static double Distance(Point2D Point1, Point2D Point2)
    {
        return Math.Sqrt(Math.Pow(Point1.X - Point2.X, 2) + Math.Pow(Point1.Y - Point2.Y, 2));
    }

    public static RelPos2D RelativePosition(Point2D Of, Point2D To)
    {
        int xRel = Of.X < To.X ? -1 : Of.X > To.X ? 1 : 0;
        int yRel = Of.Y < To.Y ? -1 : Of.Y > To.Y ? 1 : 0;

        switch (xRel)
        {
            case -1:
                switch (yRel)
                {
                    case -1: return RelPos2D.ll;
                    case 0: return RelPos2D.le;
                    case 1: return RelPos2D.lg;
                }
                break;
            case 0:
                switch (yRel)
                {
                    case -1: return RelPos2D.el;
                    case 0: return RelPos2D.ee;
                    case 1: return RelPos2D.eg;
                }
                break;
            case 1:
                switch (yRel)
                {
                    case -1: return RelPos2D.gl;
                    case 0: return RelPos2D.ge;
                    case 1: return RelPos2D.gg;
                }
                break;
        }

        return RelPos2D.ee; // never reached
    }

    public static double TriangleArea(Point2D Point1, Point2D Point2, Point2D Point3)
    {
        return 1 / 2d *
            (
                (Point1.X - Point3.X) * (Point2.Y - Point1.Y) -
                (Point1.X - Point2.X) * (Point3.Y - Point1.Y)
            );
    }

    public static bool TriangleContainsPoint(Point2D Point1, Point2D Point2, Point2D Point3, Point2D Target)
    {
        var s = Point1.Y * Point3.X - Point1.X * Point3.Y + (Point3.Y - Point1.Y) * Target.X + (Point1.X - Point3.X) * Target.Y;
        var t = Point1.X * Point2.Y - Point1.Y * Point2.X + (Point1.Y - Point2.Y) * Target.X + (Point2.X - Point1.X) * Target.Y;

        if ((s < 0) != (t < 0))
            return false;

        var area = TriangleArea(Point1, Point2, Point3);
        var sign = area < 0 ? -1 : 1;
        s *= sign;
        t *= sign;
        area *= sign;

        return s > 0 && t > 0 && (s + t) < 2 * area;
    }
}


class ProblemSolver
{
    private static RelPos2D[] AllPositions = new RelPos2D[] 
    {
        RelPos2D.ee,
        RelPos2D.eg,
        RelPos2D.el,
        RelPos2D.ge,
        RelPos2D.gg,
        RelPos2D.gl,
        RelPos2D.le,
        RelPos2D.lg,
        RelPos2D.ll,
    };

    private static RelPos2D[] NoPositions = new RelPos2D[0];

    private static RelPos2D[] ValidPositions(RelPos2D Pos1, RelPos2D Pos2)
    {
        if (Pos1 == RelPos2D.ee || Pos2 == RelPos2D.ee)
            return AllPositions;

        switch (Pos1)
        {
            case RelPos2D.ll:
                switch (Pos2)
                {
                    case RelPos2D.ll:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gg };
                    case RelPos2D.le:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gg, RelPos2D.ge };
                    case RelPos2D.lg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gg, RelPos2D.ge, RelPos2D.gl };
                    case RelPos2D.eg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gg, RelPos2D.ge, RelPos2D.gl, RelPos2D.el };
                    case RelPos2D.gg:
                        return AllPositions;
                    case RelPos2D.ge:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.le, RelPos2D.lg, RelPos2D.eg, RelPos2D.gg };
                    case RelPos2D.gl:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.lg, RelPos2D.eg, RelPos2D.gg };
                    case RelPos2D.el:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.eg, RelPos2D.gg };
                }
                break;
            case RelPos2D.le:
                switch (Pos2)
                {
                    case RelPos2D.ll:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gg, RelPos2D.ge };
                    case RelPos2D.le:
                        return NoPositions;
                    case RelPos2D.lg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ge, RelPos2D.gl };
                    case RelPos2D.eg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ge, RelPos2D.gl, RelPos2D.el };
                    case RelPos2D.gg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ge, RelPos2D.gl, RelPos2D.el, RelPos2D.ll };
                    case RelPos2D.ge:
                        return AllPositions.Except(new RelPos2D[] { Pos1, Pos2 }).ToArray();
                    case RelPos2D.gl:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.lg, RelPos2D.eg, RelPos2D.gg };
                    case RelPos2D.el:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.eg, RelPos2D.gg };
                }
                break;
            case RelPos2D.lg:
                switch (Pos2)
                {
                    case RelPos2D.ll:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gg, RelPos2D.ge, RelPos2D.gl };
                    case RelPos2D.le:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ge, RelPos2D.gl };
                    case RelPos2D.lg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gl};
                    case RelPos2D.eg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gl, RelPos2D.el };
                    case RelPos2D.gg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gl, RelPos2D.el, RelPos2D.ll };
                    case RelPos2D.ge:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gl, RelPos2D.el, RelPos2D.ll, RelPos2D.le };
                    case RelPos2D.gl:
                        return AllPositions;
                    case RelPos2D.el:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.eg, RelPos2D.gg, RelPos2D.ge, RelPos2D.gl };
                }
                break;
            case RelPos2D.eg:
                switch (Pos2)
                {
                    case RelPos2D.ll:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gg, RelPos2D.ge, RelPos2D.gl, RelPos2D.el };
                    case RelPos2D.le:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ge, RelPos2D.gl, RelPos2D.el };
                    case RelPos2D.lg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gl, RelPos2D.el };
                    case RelPos2D.eg:
                        return NoPositions;
                    case RelPos2D.gg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.el, RelPos2D.ll };
                    case RelPos2D.ge:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.el, RelPos2D.ll, RelPos2D.le };
                    case RelPos2D.gl:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.el, RelPos2D.ll, RelPos2D.le, RelPos2D.lg };
                    case RelPos2D.el:
                        return AllPositions.Except(new RelPos2D[] { Pos1, Pos2}).ToArray();
                }
                break;
            case RelPos2D.gg:
                switch (Pos2)
                {
                    case RelPos2D.ll:
                        return AllPositions;
                    case RelPos2D.le:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ge, RelPos2D.gl, RelPos2D.el, RelPos2D.ll };
                    case RelPos2D.lg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gl, RelPos2D.el, RelPos2D.ll };
                    case RelPos2D.eg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.el, RelPos2D.ll };
                    case RelPos2D.gg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ll };
                    case RelPos2D.ge:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ll, RelPos2D.le};
                    case RelPos2D.gl:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ll, RelPos2D.le, RelPos2D.lg };
                    case RelPos2D.el:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ll, RelPos2D.le, RelPos2D.lg, RelPos2D.eg };
                }
                break;
            case RelPos2D.ge:
                switch (Pos2)
                {
                    case RelPos2D.ll:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.le, RelPos2D.lg, RelPos2D.eg, RelPos2D.gg };
                    case RelPos2D.le:
                        return AllPositions.Except(new RelPos2D[] { Pos1, Pos2 }).ToArray();
                    case RelPos2D.lg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.gl, RelPos2D.el, RelPos2D.ll, RelPos2D.le };
                    case RelPos2D.eg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.el, RelPos2D.ll, RelPos2D.le };
                    case RelPos2D.gg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ll, RelPos2D.le };
                    case RelPos2D.ge:
                        return NoPositions;
                    case RelPos2D.gl:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.le, RelPos2D.lg };
                    case RelPos2D.el:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.le, RelPos2D.lg, RelPos2D.eg };
                }
                break;
            case RelPos2D.gl:
                switch (Pos2)
                {
                    case RelPos2D.ll:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.lg, RelPos2D.eg, RelPos2D.gg };
                    case RelPos2D.le:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.lg, RelPos2D.eg, RelPos2D.gg, RelPos2D.ge };
                    case RelPos2D.lg:
                        return AllPositions;
                    case RelPos2D.eg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.el, RelPos2D.ll, RelPos2D.le, RelPos2D.lg };
                    case RelPos2D.gg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ll, RelPos2D.le, RelPos2D.lg };
                    case RelPos2D.ge:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.le, RelPos2D.lg};
                    case RelPos2D.gl:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.lg };
                    case RelPos2D.el:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.lg, RelPos2D.eg };
                }
                break;
            case RelPos2D.el:
                switch (Pos2)
                {
                    case RelPos2D.ll:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.eg, RelPos2D.gg };
                    case RelPos2D.le:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.eg, RelPos2D.gg, RelPos2D.ge };
                    case RelPos2D.lg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.eg, RelPos2D.gg, RelPos2D.ge, RelPos2D.gl };
                    case RelPos2D.eg:
                        return AllPositions.Except(new RelPos2D[] { Pos1, Pos2 }).ToArray();
                    case RelPos2D.gg:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.ll, RelPos2D.le, RelPos2D.lg, RelPos2D.eg };
                    case RelPos2D.ge:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.le, RelPos2D.lg, RelPos2D.eg };
                    case RelPos2D.gl:
                        return new RelPos2D[] { RelPos2D.ee, RelPos2D.lg, RelPos2D.eg };
                    case RelPos2D.el:
                        return NoPositions;
                }
                break;
        }
        return NoPositions;
    }

    private static bool UncertainSolution(RelPos2D Pos1, RelPos2D Pos2, RelPos2D Pos3)
    {
        RelPos2D[] array = new RelPos2D[] { Pos1, Pos2, Pos3 };

        return 
            (array.Contains(RelPos2D.ll) && array.Contains(RelPos2D.gg)) ||
            (array.Contains(RelPos2D.lg) && array.Contains(RelPos2D.gl));
    }

    public Tuple<Point2D, Point2D, Point2D> Solve(Point2D Target, params Point2D[] Points)
    {
        Dictionary<Point2D, double> distanceToTarget = new Dictionary<Point2D, double>();
        Dictionary<Point2D, RelPos2D> relativePosition = new Dictionary<Point2D,RelPos2D>();
        List<int> visited = new List<int>();

        Dictionary<RelPos2D, int> countPerPosition = new Dictionary<RelPos2D, int>()
        {
           {RelPos2D.ee,0},
           {RelPos2D.eg,0},
           {RelPos2D.el,0},
           {RelPos2D.ge,0},
           {RelPos2D.gg,0},
           {RelPos2D.gl,0},
           {RelPos2D.le,0},
           {RelPos2D.lg,0},
           {RelPos2D.ll,0}
        };

        foreach (var point in Points)
        {
            distanceToTarget.Add(point, Tools2D.Distance(point, Target));
            RelPos2D position = Tools2D.RelativePosition(point, Target);
            relativePosition.Add(point, position);
            countPerPosition[position]++;
        }

        //check countPerPosition to see if there are solutions
        int pointsCount = Points.Length;
        bool noSolutions = false;
        foreach (var key in countPerPosition.Keys)
        {
            if (countPerPosition[key] == pointsCount)
            {
                noSolutions = true;
                break;
            }
        }

        noSolutions = noSolutions ||
                countPerPosition[RelPos2D.ll] + countPerPosition[RelPos2D.le] + countPerPosition[RelPos2D.lg] == pointsCount ||
                countPerPosition[RelPos2D.lg] + countPerPosition[RelPos2D.eg] + countPerPosition[RelPos2D.gg] == pointsCount ||
                countPerPosition[RelPos2D.gg] + countPerPosition[RelPos2D.ge] + countPerPosition[RelPos2D.gl] == pointsCount ||
                countPerPosition[RelPos2D.ll] + countPerPosition[RelPos2D.el] + countPerPosition[RelPos2D.gl] == pointsCount;

        if (noSolutions)
            throw new Exception("No solutions.");

        var orderedPoints = Points.OrderBy(point => distanceToTarget[point]);
        bool found = false;

        Point2D 
            Point1 = null,
            Point2 = null,
            Point3 = null;

        RelPos2D PosPoint1,
                 PosPoint2,
                 PosPoint3;

        foreach (var point1 in orderedPoints)
        {
            Point1 = point1;
            PosPoint1 = relativePosition[Point1];

            var point2Candidates = orderedPoints.Where(p => p != Point1)
                                                .OrderBy(p => distanceToTarget[p]);

            //this should not happen because we know that we have at least one solution
            if (point2Candidates.Count() == 0)
                continue;

            foreach (var point2 in point2Candidates)
            {
                Point2 = point2;
                PosPoint2 = relativePosition[Point2];

                var point3ValidPositions = ValidPositions(PosPoint1, PosPoint2);

                var point3Candidates = orderedPoints.Where(p => p != Point1 && p != Point2 && point3ValidPositions.Contains(relativePosition[p]))
                                                    .OrderBy(p => distanceToTarget[p]);

                if (point3Candidates.Count() == 0)
                    continue;

                foreach (var point3 in point3Candidates)
                {
                    Point3 = point3;
                    PosPoint3 = relativePosition[Point3];

                    //check if already visited
                    //hash subject to conflicts
                    var hash = Point1.GetHashCode() *
                               Point2.GetHashCode() *
                               Point3.GetHashCode();

                    if (visited.Contains(hash))
                        continue;

                    if (UncertainSolution(PosPoint1, PosPoint2, PosPoint3))
                    {
                        found = Tools2D.TriangleContainsPoint(Point1, Point2, Point3, Target);
                    }
                    else
                    {
                        found = true;
                    }

                    if (found)
                        break;

                    visited.Add(hash);
                }

                if (found)
                    break;
            }

            if (found)
                break;
        }

        if (found)
            return new Tuple<Point2D, Point2D, Point2D>(Point1, Point2, Point3);

        throw new Exception("No solutions.");
    }
}

class Program
{
    static void Main(string[] args)
    {
        ProblemSolver ps = new ProblemSolver();
        Random r = new Random();

        List<Point2D> points = new List<Point2D>();
        Point2D target = new Point2D()
        {
            //X = r.NextDouble() * 10,
            //Y = r.NextDouble() * 10
            X = r.Next(11),
            Y = r.Next(11)
        };

        for (int i = 0; i < 10; i++)
            points.Add(new Point2D()
            {
                //X = r.NextDouble() * 10,
                //Y = r.NextDouble() * 10
                X = r.Next(11),
                Y = r.Next(11)
            });
        Console.WriteLine("Target: {0}X: {1}{0}Y: {2}{0}", Environment.NewLine, target.X, target.Y);
        Stopwatch sw = new Stopwatch();
        sw.Start();
        try
        {
            var solution = ps.Solve(target, points.ToArray());
            Console.WriteLine("Solution: {0}X1: {1}{0}Y1: {2}{0}X2: {3}{0}Y2: {4}{0}X3: {5}{0}Y3: {6}{0}",
                Environment.NewLine,
                solution.Item1.X,
                solution.Item1.Y,
                solution.Item2.X,
                solution.Item2.Y,
                solution.Item3.X,
                solution.Item3.Y
            );
        }
        catch (Exception ex)
        {
            Console.WriteLine(ex.Message);
        }

        sw.Stop();
        Console.WriteLine("Solved in {0} ms", sw.ElapsedMilliseconds);
        Console.ReadLine();
    }