How to find the triangle
New Answer
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.
Original Answer
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.
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}}]
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.