get nearest points that form a triangle

Question

enter image description here

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:

enter image description here

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:

enter image description here

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.

Please clarify what exactly the condition is: Should the sum of the distances or the area of the triangle be minimized? — TaW, Jul 13 '14 at 12:53
@TaW: Hmmmm... Actually: I don't know. I'd like to get the "most equilateral triangle". But in fact: I think, that it doesn't matter which of this conditions is fullfilled. — Rico-E, Jul 13 '14 at 13:44
Algorithm for determining whether a point is inside a triangle available here on math.stackexchange: http://math.stackexchange.com/questions/51326/determining-if-an-arbitrary-point-lies-inside-a-triangle-defined-by-three-points — Pieter Geerkens, Jul 13 '14 at 16:10

score 2 · Answer 1 · edited May 23 '17 at 12:13

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

enter image description here

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 1

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

Equation 2

Which gives you the following formulas for line equation coefficients:

Equation 3

Equation 4

Equation 5

That will indeed find two nearest points. But not necessarily the three nearest points or the smallest triangle.. Not is it guaranteed to find any solution even if there is one.. — TaW, Jul 13 '14 at 12:21
But I like it anyway. Maybe it should be repeated with the next closest points until a solution is found. We also need to know if any solution is ok or how it must be optimised. — TaW, Jul 13 '14 at 12:29
I tested this method and I must admit that for 1000 runs it is slightly faster (2.25%) than my method, but I think I can improve it. — B0Andrew, Jul 15 '14 at 12:56
@B0Andrew Consider materializing your LINQs early. You have two `OrderBy` calls that are evaluated more than once not necessarily. — BartoszKP, Jul 15 '14 at 13:16

B0Andrew · Answer 2 · 2014-07-15T08:55:03.740

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();
    }

Sounds like an interesting idea. I was going to propose storing a list of both distances and angles and then picking from small distances and large enough angle differences. I also have a hunch that some trial and error is involved and the best solution can't be found on a direct way. But that is just a gut feeling.. — TaW, Jul 13 '14 at 16:17
I agree to TaW, this sounds like a really good idea. I'm looking forward to your example code. — Rico-E, Jul 13 '14 at 16:47
It also might help to start with the Delaunay Triangulation and angle sort the adjacent edges and find the nearest point to in the triangulation to the query point O(log N). It would probably reduce your search space on average if there are many query points on the same set of points. — Nuclearman, Jul 13 '14 at 20:43

Micromega · Answer 3 · 2014-07-13T21:02:46.680

0

You can use bowyer-watson algorithm and modify it to look for the red point.

edited Jul 13 '14 at 21:02

answered Jul 13 '14 at 20:54

Micromega

12,486
7
35
72

get nearest points that form a triangle

3 Answers3