Connecting the points surrounding the adjacent polygons

I'm trying to connect the points surrounding the adjacent polygons. I know how to find the adjacent polygons on a Voronoi mesh. However, I'm unable to adjacent of the adjacent polygons.

What I'm trying to do is to draw another line surrounding the existing black lines in the figure below.

-
Related, but not a duplicate: How to calculate a Delaunay mesh nerve? – MarcoB Mar 19 at 5:05

Here is a slightly different approach that has been packaged into an easy to call function. The idea is similar to MarcoB's approach.

pts = BlockRandom[SeedRandom["RunnyKine"]; RandomReal[5, {50, 2}]];(*for reproducibility*)

vor = VoronoiMesh[pts, MeshCellHighlight -> {{0, All} -> Blue, {1, All} -> Red,
{2, All} ->  Darker[Yellow, 0.1]}]
(* because I don't like the default appearance of VoronoiMesh *)


The following function outerBoundary takes the set of points that generates the Voronoi diagram, the Voronoi mesh itself, and the polygon of interest and returns a list containing four items: the immediate boundary of the polygon of interest, the outer boundary to this immediate boundary, the line connecting the sites of this outer boundary and the site points themselves for use in a Graphics.

 outerBoundary[pts : {{_, _} ..}, vor_MeshRegion, center_Polygon] :=
Module[{bd, firstBoundary, secondBoundary, region, rm, bpts,
boundarypoly, tiles, line},
tiles = MeshPrimitives[vor, 2];
boundarypoly[polys_][cen_] :=
Complement[
Pick[polys,
IntersectingQ[cen[[1]], #] & /@ polys[[All, 1]]], {cen}];
bd = boundarypoly[tiles];
firstBoundary = bd @ center;
secondBoundary =
Complement[bd /@ firstBoundary // Flatten, firstBoundary, {center}];
rm = RegionMember @ RegionUnion @ secondBoundary;
bpts = Pick[pts, rm[pts]];
line = Line[bpts[[Last @ FindShortestTour @ bpts]]];
{firstBoundary, secondBoundary, line, Point[bpts]}]


Here is an example usage:

tile = MeshPrimitives[vor, 2][[2]] (* random tile of interest *)

res = outerBoundary[pts, vor, tile];
Show[vor,
Graphics[{EdgeForm[Black], Red, center, Blue, res[[1]], Darker@Green,
res[[2]], Red, Thick, res[[3]], Black, PointSize[0.02], res[[4]]}]]


And another:

tile = MeshPrimitives[vor, 2][[5]];

res2 = outerBoundary[pts, vor, tile];
Show[vor,
Graphics[{EdgeForm[Black], Red, center, Blue, res2[[1]], Darker@Green,
res2[[2]], Red, Thick, res2[[3]], Black, PointSize[0.02], res2[[4]]}]]


To select the polygon with maximum number of edges just do:

polys = MeshPrimitives[vor,2];
maxnum = Max[Length @@@ polys]; (* maximum number of edges *)
maxEdgePolygon = First @ Pick[polys, UnitStep[Subtract[Length @@@ polys, maxnum]], 1]
res3 = outerBoundary[pts, vor, maxEdgePolygon];


Let's visualize it:

Show[vor,
Graphics[{EdgeForm[Black], Red, maxEdgePolygon, Blue, res3[[1]],
Darker@Green, res3[[2]], Red, Thick, res3[[3]], Black,
PointSize[0.02], res3[[4]]}]]


-
Thank you so much @RunnyKine. any suggestion how can I do the same but for the polygon with the maximum number of edges? Thanks – Abdullah Mar 19 at 7:14
@Abdullah. I'll update my answer soon, to add that. – RunnyKine Mar 19 at 7:26
I really appreciate your help. Thank you so much! – Abdullah Mar 19 at 7:27
@Abdullah, see my update. Hope it helps. – RunnyKine Mar 19 at 7:35
Thank you so much. This is really helpful and perfect – Abdullah Mar 19 at 14:51

Let's generate a few random points and get a mesh, then select one random cell in the mesh as the center for the following calculations:

SeedRandom[10]
pts = RandomReal[{-10, 10}, {70, 2}];
vm = VoronoiMesh[pts];
primitives = MeshPrimitives[vm, 2];
center = RandomChoice[primitives];


Here are mesh and center:

Show[
vm,
Graphics[{Red, center}]
]


Below I will reuse a handy function I posted earlier to find the nearest neighbors to a specified cell in a mesh (this was referred to as the "nerve"):

Clear[findnerve]
findnerve[ctr_] :=
Cases[
primitives,
Polygon[p__] /; (ContainsAny[p, ctr[[1]]] && p != ctr[[1]])
]


Your region in yellow in your example would be this "first-level nerve" nerve1:

nerve1 = findnerve[center];


In order to find the "second-level nerve", i.e. the outer-facing neighbors of those cells in nerve1, I find the neighbors of all cells in nerve1, then exclude those that were in nerve1 and the center cell:

nerve2 = Complement[
findnerve /@ findnerve[center] // Flatten // DeleteDuplicates,
findnerve[center], {center}
];

Show[
vm,
Graphics[{
Red, center,
Green, EdgeForm[Darker@Green], nerve1,
Gray, EdgeForm[Darker@Gray], nerve2
}]
]


Now that we know what polygons make up nerve2, we can select the points that belong to those polygons among the points pts used to construct the mesh, and calculate a line going through them:

nerve2points = Select[pts, RegionMember[RegionUnion[nerve2]]];
line2 = Line[nerve2points[[Last@FindShortestTour[nerve2points]]]];


Finally, we are ready to draw it all out:

Show[
vm,
Graphics[{
Red, center,
Green, EdgeForm[Darker@Green], nerve1,
Gray, EdgeForm[Darker@Gray], nerve2,
Thick, Red, line2,
PointSize[0.02], Black, Point@nerve2points
}]
]


-
Thank you so much for such a great elaboration – Abdullah Mar 19 at 5:16
In case I want to do this for the polygon with the highest number of edges. Do I need to set nerve1 to be equal to be the centre of that polygon? Thank You – Abdullah Mar 19 at 5:22
@Abdullah you are very welcome! To do the same with the polygon with the most edges, you should set nerve1 equal to the Polygon object itself, not its center. – MarcoB Mar 19 at 5:34
Thank you. polygon object as polygon's coordinates? Sorry for the trouble I'm just new to this – Abdullah Mar 19 at 6:04