I am using Voronoi polygonization as a support for simulating island erosion; my data is the output of VoronoiDiagram, DelaunayTriangulation, the starting point and the list of points that now are sea. I am not working with BoundedDiagram so I made sure all unbounded polygons are at sea.

Now I need to show this output like in DiagramPlot, but with sea cells colored blue, and land cells colored in brown.
Can you give me some code example or indications on where to search?

More generally, I need to have a planar graph of which I know vertex coordinates and the edges, what function should I use to choose the color for the faces of the graph?

up vote 12 down vote accepted

This is just a flow of consciousness... I think you are talking about something like this:

<< ComputationalGeometry`

data = {{4.4, 14}, {6.7, 15.25}, {6.9, 12.8}, {2.1, 11.1}, {9.5, 14.9}, 
       {13.2, 11.9}, {10.3, 12.3}, {6.8, 9.5}, {3.3, 7.7}, {0.6, 5.1}, 
       {5.3, 2.4}, {8.45, 4.7}, {11.5, 9.6}, {13.8, 7.3}, {12.9, 3.1}, {11, 1.1}};

ubd = DiagramPlot[data, LabelPoints -> False]

enter image description here

So here I think what you want to do:

Graphics[{FaceForm[Brown], EdgeForm[White], 
  Polygon[Select[vorvert[[#]] & /@ vorval[[All, 2]], 
    Cases[#, _Ray] == {} &]]}, Background -> Blue]

enter image description here

On the other hand image processing does it too:

MatrixPlot[MorphologicalComponents[ubd // Binarize] /. {1 -> Blue, 0 -> White, 
   x_Integer /; (x =!= 1) -> Brown}, Frame -> False]

enter image description here

but the payoff is that you lost information about Graphics primitives and went to rasterized images, - which may be fine if you care about visual only.

If you'd figure out how to sea-land with BoundedDiagram-type diagrams, then it maybe easier, because there is another way to make Voronoi. For example here coloring only approximately circular core of Voronoi cells:

plt = ListDensityPlot[RandomReal[1, {500, 3}], 
   InterpolationOrder -> 0];
plg = Cases[plt, Polygon[{x__}] :> x, ∞];
cplg = Select[plg, 
   EuclideanDistance[Mean[plt[[1, 1]][[#]]], 
      Mean[plt[[1, 1]]]] < .3 &];
Graphics[{
   GraphicsComplex[
   plt[[1, 1]], {FaceForm[Blue], EdgeForm[Black], Polygon[plg]}],
  GraphicsComplex[
   plt[[1, 1]], {FaceForm[Red], EdgeForm[White], Polygon[cplg]}]
  }]

enter image description here

  • 3
    I believe you could define plg = Cases[plt, Polygon[{x__}] :> x, ∞] (+1) – Mr.Wizard Sep 23 '13 at 13:08
  • @Mr.Wizard Yes, much better, updated, thanks ;) – Vitaliy Kaurov Sep 23 '13 at 18:03
  • :-) -- by the way, what is the benefit of keeping the GraphicsComplex form here? Since the points must be resolved for EuclideanDistance anyway it seems to me it would be simpler to use coordinates directly. – Mr.Wizard Sep 24 '13 at 1:12
  • You might consider using Nearest as belisarius did here. – Mr.Wizard Sep 24 '13 at 7:23
  • 2
    @Vitaliy Kaurov, did you have any chance to revisit it? :) – VividD May 13 '14 at 13:31

Well, another solution using V10 functionality:

data = {{4.4, 14}, {6.7, 15.25}, {6.9, 12.8}, {2.1, 11.1}, {9.5, 
    14.9}, {13.2, 11.9}, {10.3, 12.3}, {6.8, 9.5}, {3.3, 7.7}, {0.6, 
    5.1}, {5.3, 2.4}, {8.45, 4.7}, {11.5, 9.6}, {13.8, 7.3}, {12.9, 
    3.1}, {11, 1.1}};

The Voronoi diagram:

vor = VoronoiMesh[data]

Mathematica graphics

We get the Polygons that make up the Voronoi diagram and select the unbounded ones which will represent the sea, while the bounded ones will represent the land. The unbounded polygons must have at least one point in their RegionBounds identical to the overall RegionBounds of the Voronoi diagram.

cells = Map[MeshCoordinates[vor][[#]] &, MeshCells[vor, 2], {2}];

EDIT A much easier way to do the above is:

cells = MeshPrimitives[vor, 2]; (* I knew there was a better way *)

The region bounds of the Voronoi diagram is obtained simply as

RegionBounds[vor]

{{-2.7, 17.1}, {-2.4375, 18.7875}}

Using the new IntersectingQ we can see which polygons are unbounded/bounded:

bound = IntersectingQ[Flatten@RegionBounds[vor], Flatten@RegionBounds[#]] & /@ cells;

Now we can plot the regions:

polys = Transpose[{bound, cells}] /. {True -> Blue, False -> Brown};

Graphics[{EdgeForm[{White}], polys, Yellow, Point[data]}]

Mathematica graphics

Here is a nicely compact solution, which exploits the fact that the "sea cells" are precisely the cells that contain the points of the data set's convex hull. Another key ingredient is the use of the undocumented function Region`Mesh`MeshMemberCellIndex[] (previously presented by ilian here).

data = {{4.4, 14}, {6.7, 15.25}, {6.9, 12.8}, {2.1, 11.1}, {9.5, 14.9}, {13.2, 11.9},
        {10.3, 12.3}, {6.8, 9.5}, {3.3, 7.7}, {0.6, 5.1}, {5.3, 2.4}, {8.45, 4.7},
        {11.5, 9.6}, {13.8, 7.3}, {12.9, 3.1}, {11, 1.1}};

vor = VoronoiMesh[data];
ci = Region`Mesh`MeshMemberCellIndex[vor];
pos = Rest /@ ci[MeshCoordinates[ConvexHullMesh[data]]];
cells = MeshPrimitives[vor, 2];
Graphics[{EdgeForm[White], {{Blue, Extract[cells, pos]}, {Brown, Delete[cells, pos]}},
          Yellow, Point[data]}]

Voronoi sea and land

data = {{4.4, 14}, {6.7, 15.25}, {6.9, 12.8}, {2.1, 11.1}, {9.5, 
    14.9}, {13.2, 11.9}, {10.3, 12.3}, {6.8, 9.5}, {3.3, 7.7}, {0.6, 
    5.1}, {5.3, 2.4}, {8.45, 4.7}, {11.5, 9.6}, {13.8, 7.3}, {12.9, 
    3.1}, {11, 1.1}};

Since v10 there's HighlightMesh which allows to highlight the interior and boundary cells in an automated manner (see also this thread for some further details about its undocumented features). Let

vor = VoronoiMesh[data];

Then with

h = MeshCellIndex[vor, {2, "Interior"}]

{{2, 3}, {2, 4}, {2, 5}, {2, 12}, {2, 14}, {2, 16}}

and

b = MeshCellIndex[vor, {2, "Frontier"}]

{{2, 1}, {2, 2}, {2, 6}, {2, 7}, {2, 8}, {2, 9}, {2, 10}, {2, 11}, {2, 13}, {2, 15}}

one can do

highl = HighlightMesh[vor, {Style[h, Brown], Style[b, Blue]}]

enter image description here

Or also with the cell centroids:

pts = PropertyValue[{vor, 2}, MeshCellCentroid];
Show[highl, Graphics[{Red, PointSize[Large], Point[pts]}]]

enter image description here


Also, if the concern is only about the image, for fun:

points = ListPlot[data, PlotStyle -> {Red, PointSize[Large]}];

regBounds = RegionBoundary @ ConvexHullMesh @ data;

f[{x_, y_}] := If[RegionMember[regBounds, {x, y}], 1, -1]
c = ListContourPlot[Function[{x, y}, {x, y, f[{x, y}]}] @@@ data, 
     Mesh -> All, InterpolationOrder -> 0, Frame -> False,
     PlotRange -> RegionBounds @ vor];

plot = Show[
  c /. (Rule @@@ Transpose@{Cases[c, _RGBColor, Infinity], {Brown, Blue}}), 
  points]

enter image description here

MeshRegion[vor, MeshCellStyle -> {{2,_}->Yellow, {2,"Frontier"} -> Red}]

Mathematica graphics

MeshRegion[vor, MeshCellStyle -> {{2,_}->Yellow, {2,"Interior"} -> Red}]

Mathematica graphics

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.