# Quirk in VertexColors interpolation when displaying Polygon

When using VertexColors, the docs state:

The interior of the polygon is colored by interpolating between the colors specified by VertexColors.

It however doesn't say how this interpolation is carried out. I've seen at times that small quirks show up that I wouldn't expect, and I'm wondering if someone else has a clearer perspective on exactly why and when these show up.

 aim = {Cos[2 π #], Sin[2 π #]} &;
angles[n_] := Range[0, 0.5, 0.5/(n - 1)];
colors[n_] := Join[ConstantArray[Purple, n], ConstantArray[Orange, n]];

draw[n_, a_] :=
Graphics[{
Polygon[Join[aim /@ angles[n], a aim@# & /@ Reverse[angles[n]]],
VertexColors -> colors[n]
]}]


Using this definition and n = 5, I don't see anything wrong when a = 0.415, but when a = 0.414, it looks like it's interpolating incorrectly, and when n is set larger, two such quirks appear which seem to always be present. This can be seen in the three calls:

 GraphicsColumn[{draw[5, 0.415], draw[5, 0.414], draw[100, 0.415]}]


So the question is, why does this happen and is there a general way to avoid it, or test for when it's expected to happen?

I know that for this particular example I can avoid it simply by splitting my single polygon into many smaller 4 point polygons, I'm not asking for code that just draws a shaded arc. I'm interested in the underlying problem.

-
did you check if all the polygons are oriented counter clockwise? Perhaps that's an issue. – user21 Dec 12 '12 at 12:37
@ruebenko There is only one polygon. It happens to be oriented counter clockwise, but the same behavior is present independent of orientation. – jVincent Dec 12 '12 at 12:46
My guess is that the best general way to deal with this is to (1) triangulate your large polygon (2) assign vertex colors to the newly introduced vertices (could be tricky, in general) and (3) display the whole thing. – Mark McClure Dec 12 '12 at 15:31
@MarkMcClure A triangulation shouldn't introduce new vertices, unless you mean to resolve self intersections simultaneously. In such cases I would suspect it would be trivial to just do a linear interpolation at the intersection point. – jVincent Dec 12 '12 at 15:51
Depends on how you triangulate! :) I described one tool for triangulation in this answer: mathematica.stackexchange.com/questions/3646/… Of course, that's not necessary for this example but, perhaps, one would find themselves using that on another example. – Mark McClure Dec 12 '12 at 16:02

The behavior you are seeing is, that OpenGL (which is most probably used) has to break your your polygon of many vertices into triangles. What you then see are artifacts of the linear interpolation of the colors. When you assign a random color to each vertex, you can observe the triangle structure

aim = {Cos[2 π #], Sin[2 π #]} &;
angles[n_] := Range[0, 0.5, 0.5/(n - 1)];
colors[n_] :=
Join[ConstantArray[Purple, n], ConstantArray[Orange, n]];

draw[n_, a_] :=
Graphics[{
Polygon[Join[aim /@ angles[n], a aim@# & /@ Reverse[angles[n]]],
VertexColors -> (RGBColor @@@ RandomReal[{0, 1}, {2 n, 3}])]}]
Manipulate[draw[n, a], {a, 0.4, 0.5}, {n, 5, 20, 1}]


To my knowledge, it is not possible to influence the process of breaking a polygon down to triangles and therefore ensure a satisfying colorization. What you can do of course is, give triangles as polygons and therefore ensure this manually.

In your case this would be something like

points[n_, a_] := {#1, #2, #4, #3} & @@@
Partition[Flatten[Transpose[{a #, #}], 1], 4, 2] &@
Table[{Cos[φ], Sin[φ]}, {φ, 0, π, π/(n - 1)}];
Manipulate[
Graphics[Polygon[points[n, a],
VertexColors ->
ConstantArray[{Orange, Purple, Purple, Orange}, n - 1]]], {a, 0.4,
0.5}, {n, 5, 20, 1}]


-
So really the problem is that it's doing a less then optimal triangulation which actually seems to result in zero area triangles at times. Is there any documentation or information about which method is used to break it down? Perhaps someone knowledgeable in triangulation algorithms would know, Afaik Delaunay triangulation avoids these effects. – jVincent Dec 12 '12 at 13:16
Yes, I'm sure but not in Mathematica, because this is an OpenGL problem. Look for instance down this page. Looks familiar, right? If you want to know more, you could check the sources of the Mesa free OpenGL implementation. – halirutan Dec 12 '12 at 13:25

As can be seen from the images presented, the interpolation being done under the hood by VertexColors depends on a prior triangulation of the polygon, thus resulting in visible triangular bands. The approaches presented thus far have all needed to perform a preliminary triangulation; I shall now present a method that avoids this preprocessing step.

One particularly beautiful method I have encountered in the literature is the use of Floater and Hormann's mean value coordinates, which are a generalization of the classical barycentric coordinates for triangles. (There are other generalizations; if you are interested, have a look at their bibliography.) In fact, one of the applications presented in their paper is exactly this problem.

Here, now, is a Mathematica routine for performing mean value interpolation over an arbitrary polygon:

MeanValueInterpolation[poly_ /; MatrixQ[poly, NumericQ],
vals_ /; ArrayQ[vals, _, NumericQ],
x_ /; VectorQ[x, NumericQ]] /;
Dimensions[poly] == {Length[vals], Length[x]} :=
Module[{al, cq, edg, idx, k, ki, n, r, s, w},
s = # - x & /@ poly; r = Norm /@ s; idx = Position[r, w_ /; w == 0];
If[idx =!= {}, Return[Extract[vals, First[idx]]]];
cq = Boole[TrueQ[First[poly] == Last[poly]]];
edg = Partition[s/r, 2, 1, {1, 1 - 2 cq}]; n = Length[poly] - cq;
r = Take[r, n]; al = Sign[Det /@ edg]; idx = Position[al, 0];
If[idx =!= {} && Negative[Dot @@ Extract[edg, First[idx]]],
ki = Mod[idx[[1, 1]] + {0, 1}, n, 1];
Return[r[[ki]].vals[[Reverse[ki]]]/Total[r[[ki]]]]];
w = ListCorrelate[{1, 1}, al (Norm[Subtract @@ #]/Norm[Total[#]] &
/@ edg), -1]/r;
w.Take[vals, n]/Total[w]]


There are some slight differences in this implementation from the pseudocode presented in the paper; in particular, I elected to use a more stable method for computing the angle between two vectors due to Velvel Kahan.

To use this for coloring a polygonal region, we will also use RegionPlot[] along with the undocumented InPolygonQ[] for testing if a point lies entirely within a polygon.

(* for older versions; use GraphicsPolygonUtils instead in new Mathematica *)a
GraphicsMeshMeshInit[];


I'll use my own examples first, and come to the example in the OP later. Here is a simplified outline of the United States, and a list of colors corresponding to the vertices:

usa = {{-2.151, 0.858}, {-1.562, 0.862}, {-1.414, 0.769}, {-1.243, 0.774},
{-1.083, 0.846}, {-1.052, 0.787}, {-1.378, 0.653}, {-1.433, 0.547},
{-1.42, 0.461}, {-1.605, 0.521}, {-1.715, 0.461}, {-1.881, 0.563},
{-2.04, 0.578}, {-2.138, 0.629}, {-2.2, 0.738}};
cols = ColorData["ThermometerColors"] /@ Range[0, 1, 1/14];


Let's compare the mean value approach (left) and the built-in VertexColors:

GraphicsRow[{RegionPlot[InPolygonQ[usa, {x, y}],
{x, -2.2, -1.1}, {y, 0.46, 0.87},
AspectRatio -> Automatic, BoundaryStyle -> None,
ColorFunction -> (RGBColor @@ Clip[
MeanValueInterpolation[usa, List @@@ cols, {#1, #2}],
{0, 1}] &), ColorFunctionScaling -> False],
Graphics[Polygon[usa, VertexColors -> cols], Frame -> True]}]


As a second example, here is a bean-shaped polygon and a set of associated colors:

bean = {{0., 0.11}, {-0.05, 0.108}, {-0.1, 0.115}, {-0.16, 0.15}, {-0.22, 0.205},
{-0.32, 0.3}, {-0.4, 0.358}, {-0.5, 0.42}, {-0.55, 0.436}, {-0.6, 0.43},
{-0.644, 0.4}, {-0.66, 0.35}, {-0.655, 0.3}, {-0.635, 0.2}, {-0.595, 0.1},
{-0.552, 0.}, {-0.5, -0.105}, {-0.44, -0.2}, {-0.4, -0.25}, {-0.35, -0.3},
{-0.3, -0.344}, {-0.204, -0.4}, {-0.1, -0.436}, {0., -0.448},
{0.1, -0.442}, {0.23, -0.4}, {0.3, -0.35}, {0.353, -0.3}, {0.43, -0.2},
{0.477, -0.1}, {0.51, 0.}, {0.522, 0.1}, {0.52, 0.16}, {0.5, 0.24},
{0.456, 0.3}, {0.4, 0.33}, {0.36, 0.337}, {0.3, 0.32}, {0.25, 0.29},
{0.2, 0.245}, {0.15, 0.2}, {0.1, 0.16}, {0.05, 0.128}};
cols = {{0.834, 1., 0.}, {0.497, 1., 0.}, {0.172, 1., 0.}, {0., 1., 0.035},
{0., 1., 0.107}, {0., 1., 0.161}, {0., 1., 0.206}, {0., 1., 0.255},
{0., 1., 0.288}, {0., 1., 0.339}, {0., 1., 0.406}, {0., 1., 0.473},
{0., 1., 0.528}, {0., 1., 0.648}, {0., 1., 0.783}, {0., 1., 0.951},
{0., 0.834, 1.}, {0., 0.6, 1.}, {0., 0.457, 1.}, {0., 0.294, 1.},
{0., 0.137, 1.}, {0.126, 0., 1.}, {0.383, 0., 1.}, {0.612, 0., 1.},
{0.83, 0., 1.}, {1., 0., 0.887}, {1., 0., 0.714}, {1., 0., 0.571},
{1., 0., 0.332}, {1., 0., 0.133}, {1., 0.044, 0.}, {1., 0.207, 0.},
{1., 0.301, 0.}, {1., 0.43, 0.}, {1., 0.547, 0.}, {1., 0.64, 0.},
{1., 0.692, 0.}, {1., 0.745, 0.}, {1., 0.775, 0.}, {1., 0.789, 0.},
{1., 0.809, 0.}, {1., 0.853, 0.}, {1., 0.952, 0.}};


Now, color the bean in two different ways:

GraphicsRow[{RegionPlot[InPolygonQ[bean, {x, y}], {x, -0.7, 0.6}, {y, -0.5, 0.5},
AspectRatio -> Automatic, BoundaryStyle -> None,
ColorFunction -> (RGBColor @@ Clip[
MeanValueInterpolation[bean, cols, {#1, #2}], {0, 1}] &),
ColorFunctionScaling -> False],
Graphics[Polygon[bean, VertexColors -> (RGBColor @@@ cols)],
Frame -> True]}]


Finally, here is the OP's example:

aim = {Cos[2 π #], Sin[2 π #]} &; angles[n_] := Range[0, 0.5, 0.5/(n - 1)];
sector[a_, n_] := Join[aim /@ angles[n], a aim@# & /@ Reverse[angles[n]]];
colors[n_] := Join[ConstantArray[Purple, n], ConstantArray[Orange, n]];

GraphicsRow[{RegionPlot[InPolygonQ[sector[0.414, 5], {x, y}],
{x, -1, 1}, {y, 0, 1},
AspectRatio -> Automatic, BoundaryStyle -> None,
ColorFunction -> (RGBColor @@ Clip[
MeanValueInterpolation[sector[0.414, 5],
List @@@ colors[5], {#1, #2}],
{0, 1}] &),
ColorFunctionScaling -> False],
RegionPlot[InPolygonQ[sector[0.4, 30], {x, y}],
{x, -1, 1}, {y, 0, 1},
AspectRatio -> Automatic, BoundaryStyle -> None,
ColorFunction -> (RGBColor @@ Clip[
MeanValueInterpolation[sector[0.4, 30],
List @@@ colors[30], {#1, #2}],
{0, 1}] &),
ColorFunctionScaling -> False]}]


In all cases, the mean value coloring did a much better job than the built-in VertexColors. Of course, you no longer have a simple polygon, since RegionPlot[] has broken your polygon up into tiny triangles, but if the coloring is needed, this might be a small price to pay.

-

Here's a triangulation approach that might work easily and more generally. Let's start with on of your problematic examples:

drawing = draw[5,0.414]


Now, let's triangulate, after getting rid of 10^-17 type expressions.

GraphicsMeshMeshInit[];
pts = Chop[drawing[[1, 1, 1]]];
tri = PolygonTriangulate[pts]

(* Out: {{6, 4, 5}, {4, 2, 3}, {1, 9, 10}, {8, 2, 4}, {9, 2, 8},
{2, 9, 1}, {4, 7, 8}, {4, 6, 7}} *)


Finally, we'll assemble it into a GraphicsComplex.

colors = VertexColors /. drawing[[1, 1, 2]];
newcolors = colors[[#]] & /@ tri;
Graphics[GraphicsComplex[pts,
Polygon[tri, VertexColors -> newcolors]]]


We can see the triangulation using EdgeForm, which might be useful to understand the colorization.

 Graphics[GraphicsComplex[pts,
{EdgeForm[Black], Polygon[tri, VertexColors -> newcolors]}]]


Which explains why it doesn't work so well here with draw[100,0.414]:

On the other hand we could use the following triangulation of your draw[100,0.414] polygon (obtained from the triangle program as described in this answer):

We could then color the newly introduced vertices using linear interpolation between orange and purple. The result is:

Which is not too bad! I'd be glad to share more specific code, if you get the triangle command working, as described in the globe answer.

-
@jVincent Just wondering if you had noticed the edits to this answer? – Mark McClure Dec 12 '12 at 19:22
It would be nice to see some code pieces, you can always just reference your other answer with respect to interfacing triangle. – jVincent Dec 13 '12 at 9:41
@jVincent Yes, it would be. I was working on a bunch of stuff yesterday and figured I'd only prioritize this if you had triangle working. I'll probably clean this up at some point in the next couple of days but don't have the time immediately. I am curious to try it on more complicated examples, anyway. – Mark McClure Dec 13 '12 at 15:04

As stated in my answer, the problem for the arc was solved with just creating polygons such that they ended up without defects similar to how halirutan did it, which is a method to avoid these problems in specific cases. But in the general case this doesn't work out, so I also found a solution by triangulating the polygon using the Delaunay triangulation. This seems to work, and my intuition tells me it should work for all non-self-intersecting polygons (* My intuition fails me, see update):

 Needs["ComputationalGeometry"]

trisFromDelaunay[o_,l_]:=Sequence@@({o,Sequence@@#}&/@Partition[l,2,1])
triangulate[points_]:=trisFromDelaunay@@@DelaunayTriangulation[points]
filter[indices_]:=Select[indices,Signature[#]==1&]

polygonColoredTris[points_,VertexColors->colors_]:=
Polygon[points[[#]]&/@#,VertexColors->colors[[#]]]&/@filter[triangulate[points]]


This seems to work nicely for my current uses.

n = 50;
l = 0.3; u = 1.9;
spiral = Join[#, Reverse[0.6 #]] &@(# aim[#] & /@ Range[l, u, (u - l)/(n - 1)]);
spiralColors = Blend[{Purple, Orange}, #] & /@ Range[0, 1, 1/(n - 1)];

GraphicsRow@Table[
Graphics[{
testing[spiral, VertexColors -> colors[n]]
}],
{testing, {Polygon, polygonColoredTris}}]
`

Update Whoops, turns out this also only works under special conditions. Specifically, you aren't guaranteed that the original edges are actually present in the triangulation; I got by because I was looking at cases where they did. So I'm back to the drawing board.

-