# How to find point within polygons and connect them

I have the following figure:

Here, I have the following values of 25 image corner-points:

a = {{528.5, 563.5}, {282.5, 726.5}, {283.5, 722.5}, {42.5,
607.5}, {264.5, 394.5}, {269.5, 322.5}, {360.5, 334.5}, {366.5,
375.5}, {239.5, 324.5}, {265.5, 332.5}, {210.5, 326.5}, {366.5,
315.5}, {238.5, 394.5}, {195.5, 441.5}, {379.5, 333.5}, {219.5,
397.5}, {371.5, 386.5}, {121.5, 635.5}, {350.5, 335.5}, {350.5,
386.5}, {291.5, 320.5}, {344.5, 375.5}, {445.5, 310.5}, {342.5,
317.5}, {326.5, 390.5}}


Then, I calculated the Voronoi mesh polygons for these corner-points. The Voronoimesh polygons are:

vm= {Polygon[{{235.086, 434.547}, {123.376, 373.614}, {223.462,
360.927}}], Polygon[{{293.659, 363.962}, {304.17, 526.888}, {251.5,
449.564}, {251.5, 363.282}}],
Polygon[{{270.045, 206.5}, {281.557, 333.123}, {254.418,
322.267}, {246.7, 206.5}}],
Polygon[{{361., 523.5}, {338.104, 386.125}, {355.5, 376.636}, {361.,
384.636}}],
Polygon[{{117.5, 830.5}, {190.486, 701.371}, {492.956,
776.989}, {528.413, 830.5}}],
Polygon[{{353.98, 562.844}, {492.956, 776.989}, {190.486,
701.371}, {257.644, 576.318}, {316.267, 557.959}}],
Polygon[{{354.303, 323.032}, {354.984, 322.311}, {369.628,
326.935}, {371.047, 353.896}, {357.587, 355.865}}],
Polygon[{{428.532, 368.005}, {405.303, 364.499}, {371.047,
353.896}, {369.628, 326.935}, {405.254, 301.206}}],
Polygon[{{354.303, 323.032}, {357.587, 355.865}, {355.5,
356.7}, {315.239, 350.661}, {318.182, 339.086}}],
Polygon[{{311.529, 354.235}, {315.239, 350.661}, {355.5,
356.7}, {355.5, 376.636}, {338.104, 386.125}}],
Polygon[{{345.333, 206.5}, {354.984, 322.311}, {354.303,
323.032}, {318.182, 339.086}, {310.382, 206.5}}],
Polygon[{{-79., 830.5}, {-79., 344.139}, {-72.4138,
348.076}, {117.025, 522.68}, {7.92405, 830.5}}],
Polygon[{{399.259, 206.5}, {405.254, 301.206}, {369.628,
326.935}, {354.984, 322.311}, {345.333, 206.5}}],
Polygon[{{7.92405, 830.5}, {117.025, 522.68}, {257.644,
576.318}, {190.486, 701.371}, {117.5, 830.5}}],
Polygon[{{246.7, 206.5}, {254.418, 322.267}, {242.944,
359.556}, {227.333, 359.333}, {216.793, 206.5}}],
Polygon[{{355.5, 376.636}, {355.5, 356.7}, {357.587,
355.865}, {371.047, 353.896}, {405.303, 364.499}, {361.,
384.636}}],
Polygon[{{295.094, 362.453}, {293.659, 363.962}, {251.5,
363.282}, {242.944, 359.556}, {254.418, 322.267}, {281.557,
333.123}}],
Polygon[{{235.086, 434.547}, {223.462, 360.927}, {227.333,
359.333}, {242.944, 359.556}, {251.5, 363.282}, {251.5,
449.564}}],
Polygon[{{363.508, 551.719}, {361., 523.5}, {361.,
384.636}, {405.303, 364.499}, {428.532, 368.005}, {496.269,
433.959}}],
Polygon[{{650., 206.5}, {650., 383.526}, {496.269,
433.959}, {428.532, 368.005}, {405.254, 301.206}, {399.259,
206.5}}],
Polygon[{{528.413, 830.5}, {492.956, 776.989}, {353.98,
562.844}, {363.508, 551.719}, {496.269, 433.959}, {650.,
383.526}, {650., 830.5}}],
Polygon[{{310.382, 206.5}, {318.182, 339.086}, {315.239,
350.661}, {311.529, 354.235}, {295.094, 362.453}, {281.557,
333.123}, {270.045, 206.5}}],
Polygon[{{216.793, 206.5}, {227.333, 359.333}, {223.462,
360.927}, {123.376, 373.614}, {-72.4138, 348.076}, {-79.,
344.139}, {-79., 206.5}}],
Polygon[{{316.267, 557.959}, {257.644, 576.318}, {117.025,
522.68}, {-72.4138, 348.076}, {123.376, 373.614}, {235.086,
434.547}, {251.5, 449.564}, {304.17, 526.888}}],
Polygon[{{361., 523.5}, {363.508, 551.719}, {353.98,
562.844}, {316.267, 557.959}, {304.17, 526.888}, {293.659,
363.962}, {295.094, 362.453}, {311.529, 354.235}, {338.104,
386.125}}]}


Then, I choose 9 polygons from here. For your convenience, I am only giving the polygon values. The polygons are:

  aa = {{{293.659, 363.962}, {304.17, 526.888}, {251.5, 449.564}, {251.5,
363.282}}, {{361., 523.5}, {338.104, 386.125}, {355.5,
376.636}, {361., 384.636}}, {{353.98, 562.844}, {492.956,
776.989}, {190.486, 701.371}, {257.644, 576.318}, {316.267,
557.959}}, {{311.529, 354.235}, {315.239, 350.661}, {355.5,
356.7}, {355.5, 376.636}, {338.104, 386.125}}, {{295.094,
362.453}, {293.659, 363.962}, {251.5, 363.282}, {242.944,
359.556}, {254.418, 322.267}, {281.557, 333.123}}, {{363.508,
551.719}, {361., 523.5}, {361., 384.636}, {405.303,
364.499}, {428.532, 368.005}, {496.269, 433.959}}, {{528.413,
830.5}, {492.956, 776.989}, {353.98, 562.844}, {363.508,
551.719}, {496.269, 433.959}, {650., 383.526}, {650.,
830.5}}, {{310.382, 206.5}, {318.182, 339.086}, {315.239,
350.661}, {311.529, 354.235}, {295.094, 362.453}, {281.557,
333.123}, {270.045, 206.5}}, {{316.267, 557.959}, {257.644,
576.318}, {117.025, 522.68}, {-72.4138, 348.076}, {123.376,
373.614}, {235.086, 434.547}, {251.5, 449.564}, {304.17, 526.888}}}


Now, I want to find the corner-points from "a" which are within the polygons described in aa. Then connect them by a line. So, the final figure would be like this:

Please let me know, how to do a automatic code of that in mathematica. Thanks.

• mathematica.stackexchange.com/q/9405/2079 Mar 4, 2016 at 20:03
• Which one? Graphics /@ Polygon /@ Permutations@RandomInteger[{1, 30}, {6, 2}] Mar 4, 2016 at 21:14
• Hi, thanks for the comment. I don't get it. In the link there are lots of codes there. Should I use " inpolyQ" function that describe in the code? Mar 4, 2016 at 21:47
• You already asked eight questions, never accepted an answer and voted only once. Please be a nice site citizen. Read the following comment Mar 4, 2016 at 21:50
• Hi, Sorry. I didn't know that I have to do that. I always give thanks comment. I voted them all now. and Also accepted the correct answers. Thanks a lot. Mar 4, 2016 at 22:09

Using the FindAdjPoly function presented here we can find those cells adjacent to our one of interest. Then we can look back and find the points which were the seeds for those cells (in order thankfully) and simply connect them. I've modified the output of FindAdjPoly slightly to keep the rotational ordering. There also may be a more elegant way to retrieve the seed points.

    pts = {{528.5, 563.5}, {282.5, 726.5}, {283.5, 722.5}, {42.5,
607.5}, {264.5, 394.5}, {269.5, 322.5}, {360.5, 334.5}, {366.5,
375.5}, {239.5, 324.5}, {265.5, 332.5}, {210.5, 326.5}, {366.5,
315.5}, {238.5, 394.5}, {195.5, 441.5}, {379.5, 333.5}, {219.5,
397.5}, {371.5, 386.5}, {121.5, 635.5}, {350.5, 335.5}, {350.5,
386.5}, {291.5, 320.5}, {344.5, 375.5}, {445.5, 310.5}, {342.5,
317.5}, {326.5, 390.5}};

Block[{pt, mpt, regs, lines, poly, loc, all},
mpt = MeshCoordinates[Vmesh];
(*Get the mesh polygons*)regs = MeshCells[Vmesh, 2];
pt = pts[[ptIndex]];
(*Select the polygon containing pt*)
poly = SelectFirst[(regs /. Polygon[a__] :> {a, Polygon@mpt[[a]]}),
RegionMember[Last@#, pt] &];
(*Find the index of the polygon in mesh*)
loc = MeshCellIndex[Vmesh, {Polygon[First@poly]}][[1, 2]];
(*Get the edges of the polygon*)

lines = Partition[First@poly, 2, 1, 1];
(*Find polygons that share the above edges*)

all = DeleteDuplicates@Flatten[With[{ed = #},
Position[regs /. Polygon[a__] :> (MemberQ[a, #] & /@ ed),
{True, True}]] & /@ lines];

{loc, DeleteCases[all, loc]}
]


Then

    iCell = 25;

Vmesh = VoronoiMesh[pts];
{iPoly, polys} = FindAdjPoly[iCell, {pts, Vmesh}];

mpt = MeshCoordinates[Vmesh];
regs = MeshCells[Vmesh, 2];
cpts = Table[
SelectFirst[pts, RegionMember[(regs[[i]] /.
Polygon[a__] :> {a, Polygon@mpt[[a]]})[[2]], #] &],
{i, polys}];

Show[
HighlightMesh[Vmesh, {{{2, iPoly}},
Transpose[{ConstantArray[2, Length[polys]], polys}]}],
Graphics[{
{Blue, Point[pts]},
{Black, Point[cpts], JoinedCurve[Line[cpts], CurveClosed -> True]}
}]
]


Works for any Voronoi mesh:

    pts = RandomReal[{0, 1}, {100, 2}];
iCell = 50;
...


• Hi, Thanks for the code. But, when I change the points values to 50 instead of 25, it gives error. I changed all 25 to 50 in the code. But, it's still giving error. What should I do? Mar 4, 2016 at 23:43
• You don't have 50 seed points or polygons so of course it can't highlight and ring the 50th one; you only have 25. Curiously the 25th point seemed to be a lucky case where the 25 element in pts corresponded to the 25th polygon in Vmesh, this is in general not the case and I have amended the code above such that now changing the value of iCell will highlight and ring the corresponding polygon correctly. Mar 5, 2016 at 12:06
• Hi, thanks for the code. Now, I want to change the icell value to the value that hold the maximum adjacent point. I have that index value. But, when I putt that in the code, it gives an error. I found it it's because there is a {} in the value. Here is the code for founding the polygon that held the maximum adjacent points: pv = MeshPrimitives[Vmesh, 2] numv = Length /@ Level[pv, {2}] sv = Max[numv] nv = Count[numv, sv] index2 = Level[Position[numv, sv], {2}] Now, this index2 is giving {50}. When I put this, it gives an error. How to solve this? Mar 6, 2016 at 0:53
• If you only care about the cell with the highest number of adjacent polygons you could have put this in you question statement. Anyhow, you should have good for for doing that from this answer to your previous question: mathematica.stackexchange.com/a/105216/6588 If your problem is that you have iCell={50} instead of iCell=50 there are many ways of removing the List head, as described here: mathematica.stackexchange.com/questions/109099/… mathematica.stackexchange.com/questions/19277/… Mar 6, 2016 at 10:29
• iCell = First@index2 iCell = Sequence @@ index2 iCell = index2[[1]] iCell = Delete[index2, 0] Mar 6, 2016 at 10:31