# How to find the triangle

I have some points from a triangle,but there are some unexpected point still.How to remove those singular point and find the most probable triangle?

SeedRandom[8]
triPoint =
RandomPoint[
Triangle[{{13/20, 17/100}, {13/50, 37/100}, {3/5, 51/100}}], 200];
singular =
RandomPoint[
RegionDifference[
Polygon[{{9/50, 33/50}, {4/25, 1/10}, {41/50, 3/100}, {43/50, 57/
100}}], Triangle[{{13/20, 17/100}, {13/50, 37/100}, {3/5, 51/
100}}]], 10];


Suppose I have these points like following:

myPoints = Join[singular, triPoint];


We can show the point:

Graphics[Point[myPoints]]


Every real number is between 0 and 1.How to find the finally triangle:Triangle[{{13/20,17/100},{13/50,37/100},{3/5,51/100}}]?

## Update

This is a problem to simulate how to get a ternary phase.Well,I think use graphy theory maybe can help.Conneting all of points by this post.

ConnectSeparateGraphCostMin[graph_] :=
Module[{rule, cc, minWeightInCom, allMinEdge, spanTreeEdge,
completeGraph},
rule = Dispatch[
cc = ConnectedComponents[graph];
minWeightInCom =
First@MinimalBy[#, Last] & /@
Map[{#, EuclideanDistance @@ (# /. rule)} &,
Tuples /@ Subsets[cc, {2}], {2}];
allMinEdge = UndirectedEdge @@@ First /@ minWeightInCom;
spanTreeEdge =
EdgeList[
FindSpanningTree[
completeGraph =
CompleteGraph[Length@cc, EdgeWeight -> Last /@ minWeightInCom]]];
allMinEdge[[EdgeIndex[completeGraph, #] & /@ spanTreeEdge]]]]

connectGraph =
ConnectSeparateGraphCostMin[NearestNeighborGraph[myPoints]]


But I'm lost in how to do next step.The GraphPeriphery strike to me,but it's seem not work.

• Is it generally the case in your data that the points in the triangle are much denser than those outside? If so, you can consider doing a Voronoi tessellation and deleting points with large Voronoi area. – yohbs Dec 6 '16 at 14:27
• Are the number of points for the triangle and random selection always 200 and 10, respectively? I ask because it's not clear if you're asking for a specific solution or a general solution. Might you have 200 and 200 points, respectively? (And I'm not aware of the definition of the term "most probable triangle".) – JimB Dec 6 '16 at 19:11
• I have voted to close the question because (1) it is not obvious as to what knowledge one has to provide a solution besides the collection of data points (Does the solver know the random selection process? Does the user know the number of points belonging to the triangle?) and (2) it is not clear how one would generalize the solution with the specific problem being too narrow to help anyone else. – JimB Dec 6 '16 at 22:57
• @yohbs you mean this? – yode Dec 7 '16 at 1:44
• @JimBaldwin Done! :) – yode Dec 7 '16 at 3:44

OK, here is an almost failsafe way. If a SmoothKernelDistribution gives you a good estimate of the region where your dense triangle points lie, then this way will find the correct answer. I'm showing the behavior on a random set of points without fixing the Seed.

We use the kernel distribution to decide which points lie inside the triangle.

selectInnerPoints[pts_] := Module[{sk},
sk = SmoothKernelDistribution[pts];
Select[pts, PDF[sk, #] > 5 &]
];

points = Join[triPoint, singular];
ListPlot[{points, selectInnerPoints[points]}]


What we are going to do is simple. We find the smallest triangle that does not contain any of the inside points. Sounds simple and is simple. First, we need a function to check whether a point is inside a given triangle. Bluntly stolen from Eric and even written with the same variable names:

insideQ[{p1_, p2_, p3_}, v_] :=
Module[{v0 = p1, v1 = p2 - p1, v2 = p3 - p1, a, b},
a = (Det[{v, v2}] - Det[{v0, v2}])/Det[{v1, v2}];
b = -(Det[{v, v1}] - Det[{v0, v1}])/Det[{v1, v2}];
a > 0 && b > 0 && a + b < 1
]


Next, we just need the area of a triangle, because this is what we want to minimize. You know this from school, but in case you don't, Eric has the answer again

area[{p1_, p2_, p3_}] := 1/2 Abs[Det[{p2 - p1, p3 - p1}]]


The next would be a one-liner if it wasn't so tedious to write the variables down. We minimize the area of a starting triangle under the condition that no inner point is outside. We start with a triangle large enough to contain all points which can easily be constructed from the bounding box of all your points.

findTriangle[pts_] := Module[{xmin, xmax, ymin, ymax, innerPoints},
innerPoints = selectInnerPoints[pts];
{{xmin, xmax}, {ymin, ymax}} = MinMax /@ Transpose[pts];

FindMinimum[{area[{{p1x, p1y}, {p2x, p2y}, {p3x, p3y}}],
And @@ (insideQ[{{p1x, p1y}, {p2x, p2y}, {p3x, p3y}}, #] & /@
innerPoints)
},
{{p1x, xmin}, {p1y, ymin}, {p2x, xmax + ymax},
{p2y, ymin}, {p3x, xmin}, {p3y, xmax + ymax}}
]
]


Let's try it.

Show[
ListPlot[{triPoint, singular}],
Graphics[{EdgeForm[Black], FaceForm[Opacity[.2, Blue]],
Polygon[{{p1x, p1y}, {p2x, p2y}, {p3x, p3y}} /. res]}],
Frame -> True,
FrameTicks -> False,
Axes -> False
]


Worked for several different random point sets. I hope I deserve the accept with this.

Since you are speaking of the

the most probable triangle

Why not trying a kernel density estimator which gives you an estimation of the density in an instant? Without adjusting any parameters, let's try a quick hack. I create a SmoothKernelDistribution from your points, and as you can see, it gives a pretty good estimate where your point density is high

sk = SmoothKernelDistribution[Join[singular, triPoint]];
With[{dens = DensityPlot[
Log@(.01 + PDF[sk, {x, y}]), {x, 0, 1}, {y, 0, 1},
MaxRecursion -> 2, PlotPoints -> 50],
lp = ListPlot[Join[singular, triPoint]]
},
Manipulate[
Show[
dens,
ContourPlot[p == PDF[sk, {x, y}], {x, 0, 1}, {y, 0, 1},
ContourStyle -> Red],
lp
],
{p, 1, 10}
]
]


Now, you can use this to find 3 points that are most far apart. You can use the simple Euclidean distance between the 3 points for that.

dist[{{x1_?NumericQ,y1_},{x2_,y2_},{x3_,y3_}}]:=
(x1-x2)^2+(x1-x3)^2+(x2-x3)^2+(y1-y2)^2+(y1-y3)^2+(y2-y3)^2


And now you try to find the 3 points that maximize this distance function

prob = 3;
res = Quiet@
Last@NMaximize[{dist[{{x1, y1}, {x2, y2}, {x3, y3}}],
PDF[sk, #] > prob & /@ {{x1, y1}, {x2, y2}, {x3, y3}}}, {x1, y1,
x2, y2, x3, y3}];


Plotting the result

Show[
DensityPlot[
Log@(.01 + PDF[sk, {x, y}]), {x, 0.2, .8}, {y, 0.1, .6},
MaxRecursion -> 2, PlotPoints -> 50],
ListPlot[{triPoint, singular}],
Graphics[{EdgeForm[Red], FaceForm[],
Polygon[{{x1, y1}, {x2, y2}, {x3, y3}}] /. res}]
]


Not perfect, but not bad either.

• I think that's a good way to go (SmoothKernelDistribution). One could sort the sample points by the estimated density and exclude the 10 points with the lowest estimated densities which would be an automatic way to get the contour you find using Manipulate. But that's only if one knows that there are exactly 200 triangle points and exactly 10 non-triangle points. – JimB Dec 7 '16 at 17:46
• Yes, my answer was rather an idea that one can be used as a start. – halirutan Dec 7 '16 at 17:57
• Thanks for so constructive answer. :) – yode Dec 7 '16 at 18:02
• @halirutan thanks for this very instructive answer. Learned a lot. :) – ubpdqn Dec 9 '16 at 1:02
• Thanks. I was just throwing an idea in my initial post and did not expect the accept. I felt I had to add a bit more. – halirutan Dec 9 '16 at 1:04

Here is a first attempt. (i) collecting points close too each other (ii) determining the convex hull (iii) choosing triplet of boundary points that maximizes enclosure of points

findtg[pts_, thr_] :=
Module[{di = DistanceMatrix[pts], pos, ctg, cnv, cand, f, crit, max},
pos = Position[di, _?(# < thr &)];
ctg = pts[[Union[
Join @@ DeleteCases[
Union[pos, SameTest -> (#1 == Sort@#2 &)], {x_, x_}]]]];
cnv = ConvexHullMesh[ctg];
cand = (Triangle /@
Subsets[Join @@ (List @@@ MeshPrimitives[cnv, 0]), {3}]);
f[x_] := N[Total@Boole[RegionMember[x, ctg]]/Length[ctg]];
crit = f /@ cand;
max = Max[crit];
Pick[cand, # == max & /@ crit][[1]]]


Visualizing:

Manipulate[
Show[Graphics[{EdgeForm[Black], FaceForm[None],
Triangle[{{13/20, 17/100}, {13/50, 37/100}, {3/5, 51/100}}],
EdgeForm[Red], FaceForm[LightRed], findtg[myPoints, threshold]}],
ListPlot[myPoints], Frame -> True,
PlotLabel ->
Column[{Row[{"threshold:", threshold}], findtg[myPoints, threshold]
}]],
{threshold, {0.01, 0.02, 0.025, 0.03, 0.035, 0.04}}]


• Thanks,I post another short answer from your answer.If you mind that,I will delete it. :) – yode Dec 7 '16 at 1:42

This is an extended comment to highlight the need to assess the proposed algorithms with other sets of data. If we start with the same random number seed but then call @yode 's and @halirutan 's algorithm 21 times we see the following:

Yode:

Halirutan:

We see that there is generally very good "visual" agreement with the "true" triangle (outlined in red) and the estimated triangle (outlined in black).

However, there are instances in each that don't match what one would come up with visually. For the Yode algorithm there are two estimated triangles with much larger areas included. For the halirutan algorithm there are 3 instances where just a line segment (rather than a predicted triangle) is produced.

Whether or not these lack-of-fit rates are worth worrying about is a subject matter issue.

• This assess remind me should adjust the method to get that conv.Thanks. – yode Dec 7 '16 at 21:16

Almost!

## Find the shortest distance in every point

pointMinDist =
Association[
Rule[#, EuclideanDistance[#,
First@Nearest[Complement[myPoints, {#}], #]]] & /@ myPoints];


## Find all triangles which contain all points in a conv polygon

conv = ConvexHullMesh[
Keys[Select[pointMinDist, # < Mean[pointMinDist] &]]];
cnvPoint = Sequence @@@ MeshPrimitives[conv, 0];
tris = Select[Subsets[InfiniteLine @@@ MeshPrimitives[conv, 1], {3}],
And @@ RegionMember[
Triangle[Sequence @@@ RegionIntersection @@@ Subsets[#, {2}]],
cnvPoint] &];


## Select the min triangle

Graphics[{ternary =
First[MinimalBy[
Triangle[Sequence @@@ RegionIntersection @@@ Subsets[#, {2}]] & /@
tris, Area]], Blue, PointSize[.008],
Point[inside = Select[myPoints, RegionMember[ternary]]],
PointSize[.008], Red, Point[Complement[myPoints, inside]],
FaceForm[], EdgeForm[Red],
Triangle[{{13/20, 17/100}, {13/50, 37/100}, {3/5, 51/100}}]}]


ps:Actually we can add a coefficient to adjust the result,I just for the code pure to omit.

• it is your answer keep it...very nice – ubpdqn Dec 7 '16 at 1:43