Streaks while using color indexing with fill

I am trying to fill the area between two curves with a color that depends on a value at each point along one of the curves, but I am getting variations in color that should not be there. I have two (100x1) vectors: y1 and y2 which I plot against a third vector x (also 100x1) such that every point in y2 is directly above (or below or at the same location as) a point on y1. I have a fourth vector v that is the same length as x, y1, and y2, which I would like to use as a color index to fill the area between the curves defined by (x,y1) and (x,y2). (The datasets are attached below, for reference.) The command fill needs a closed loop to define the area to be filled, so I have created vectors X and Y to form the loop:

X=[x;flipud(x)];
Y=[y1;flipud(y2)];

If I type

fill(X,Y,'b')

the area between y1 and y2 will be filled in with blue. I want to fill this area with colors defined by the third variable, so I created C to serve as the color index (per the helpfile on fill):

C=[v;flipud(v)];
fill(X,Y,C)

I therefore expect the area between y1 and y2 to be filled with a color that varies horizontally (from one x value to the next), but not vertically, since each value of C corresponds to a single value in X. I get substantial vertical variation, however:

The command fill works as described in the documentation when I do not use a color index, so I believe that I am using it correctly. I have verified that the color index C has the same value at both the upper and lower points in x, since that would account for at least some variation in C between the curves at the same x position. There was no error in the color index, and it would not account for the degree of value variation in any event.

Is this a known issue while working with fill? Does anyone have suggestions as to 1) why this is happening and/or 2) how I can fix it?

[Sarah]

Here's the dataset I have been using-- I thought I had attached it to my original question.

Answer 1

[Sarah]

I solved this two ways: The 'quick and dirty' way, which involved plotting the transpose of my data so that the interpolation looked reasonable (even though the underlying problem was still present--see the discussion below for details), and using rotated text arrow annotations (because textboxes are not rotatable) in conjunction with getframe() to save the figure as it appeared before rotating. (This is fiddly to set up, since you have to manually set the position of the axis and tick labels. It's neither particularly quick nor particularly imprecise, but it's a faster-running fix than the next one.)

And the technically more-correct method, where I selected each consecutive pair of x locations on each curve (so, four points: two each on the upper and lower curves) to form a trapezoid, then filled the trapezoid with the value in the upper left corner (I could have done a simple arithmetic mean of the values or any other combination you could name, but applied the KISS principle instead). This was slower to run, but (surprisingly) faster to set up. It also has the advantage of being able to represent three variables in RGB space (e.g., concentrations in a three-part solution).

All of which will hopefully become obsolete as the older versions of MATLAB are used less and less.

[Mike Garrity]

Glad you found a solution.

The breaking it into trapezoids is an instance of the meshing approach that I mentioned in that blog post. If your trapezoids are small enough, you could probably use FaceColor = 'interp' with them.

The "faster to setup" part is probably due to the fact that you know more about the polygon than the patch object does. When it triangulates the patch, it's using a very general purpose triangulation algorithm. A more tailored algorithm, like the one you're using, can often be faster.

I'm not sure what the "slower to run" part is though. Perhaps you're creating lots of separate patch objects? If so, you could probably combine them into a single patch as I show in the cylinder example in this blog post.

Answer 2

[Star Strider]

I don’t have your data, but from my experience with fill, patch, and their friends, you may want to horizontally concatenate these rather than vertically concatenate them:

X=[x  flipud(x)];
Y=[y1 flipud(y2)];

[Sarah]

Thanks for your response! I just tried your solution-- transposed the datasets ( x=x', etc.) and plotted them, then did the concatenation (using fliplr instead of flipud due to the transposition) and fill, and there was, unfortunately, no change to the pattern of colors in the filled area.

What did work--at least in terms of getting the colors to behave--is transposing the plot itself:

   %Original:
   plot(x,y1), hold on, plot(x,y2) %correct curves
   fill(X,Y,C)                     %incorrect fill color pattern
   %New:
   plot(y1,x), hold on, plot(y2,x) %incorrect (transposed) curves
   fill(Y,X,C)                     %correct fill pattern (relative to curves)

Unfortunately, simply plotting it then changing the view to undo the transpose (using view(90,-90)) results in the original (incorrect) fill color pattern, and I really do need the (original) y-axis to be vertical.

Do you (or anyone else who might read this) have sufficient experience with fill and patch etc. to be able to provide guidance on why they preferentially interpolate along the horizontal?

[Sarah]

I found a blog post about it here:

http://blogs.mathworks.com/graphics/2016/01/14/polygon-interpolation/

and I'm going to try to understand it. I still can't believe that there isn't a way to make the vertical more dominant than the horizontal for interpolation.

[Star Strider]

Without at least a sample of your data to work with, it is difficult to experiment to see what works with it. The fliplr is correct, but for my example to work, the vectors both need to be row vectors.

Change my code to:

X=[x(:)'  flipud(x(:)')];
Y=[y1(:)' flipud(y2(:)')];

and see if that works.

[Mike Garrity]

What version of MATLAB are you using?

When I wrote that blog post about polygon interpolation, I was using the graphics system that was introduced in R2014b. The earlier one had a number of bugs in this area. You may have hit one of those. I'm guessing from the 'jet' colormap that you may be using an earlier version.

Here's my attempt to replicate your example with R2015b:

x = linspace(0,2*pi,100)';
y1 = sin(x);
y2 = 1+sin(x);

xp = [x;flipud(x)];
yp = [y1;flipud(y2)];
cp = xp;
fill(xp,yp,cp)

xp = [x;flipud(x)];
yp = [y1;flipud(y2)];
cp = yp;
fill(xp,yp,cp)

But both of those cases appear to work in R2014a:

So I expect I haven't reproduced your example correctly.

[Sarah]

@StarStrider: I tried transposing the vectors within the concatenation step as you describe before plotting and filling as well as first transposing, then concatenating (x=x'; X=[x fliplr(x)];) before plotting and filling. In both cases, the plotted (and filled) vectors were row vectors (of size 1x200). In both cases, the interpolated fill was identical (or at least qualitatively so) to the very first case (with no transposition at all).

@Mike Garrity The version I am using is indeed R2014a. Your blog post indicates that color interpolation goes all right as long as you indexing to a linear function of X and/or Y-- but my color indexing vector is not linearly related to either. In your cases showing vertical coloration (that is, the first of each of your pairs), you explicitly set cp=xp, which is about as simply and linearly related as it's possible to get!

What makes this very, very strange is that the interpolation is darn-near perfect if I plot my X values as y's and viceversa:

fill(Y,X,C)    % transposed curves

yields:

(I used the lines colormap for greater contrast/visibility)

So I'm working through your blog post (somewhat hampered by the program insisting that it has the pde toolbox and a licesnse for it, but no knowledge of createpde but I will persevere) since I cannot use the trick of using a linear relationship between C and X or the meshgrid workaround.

[Star Strider]

The only other possibility (that just now occurred to me) is that your vectors aren’t sorted by ‘x’. Concatenate your data into a single (Nx2) matrix, sort them (using the sortrows function, sorting by the column with the ‘x’ data) and see if doing that helps. I have no idea what your data are, so I can’t test this idea.

[Mike Garrity]

@Mike Garrity The version I am using is indeed R2014a. Your blog 
post indicates that color interpolation goes all right as long as 
you indexing to a linear function of X and/or Y-- but my color 
indexing vector is not linearly related to either.

Well then, as I explained in that blog post, the patch object is not going to be able to do this for you.

somewhat hampered by the program insisting that it has the pde 
toolbox and a licesnse for it ...

I was just using pdetool to illustrate the approach because it takes care of a lot of the grunt work for you. The basic idea is the same regardless of how you do it.

If the function describing your colors is not a linear combination of your X & Y coordinates, then you are going to have to subdivide your shape into enough triangles that the linear interpolation within each triangle is a reasonably accurate representation of your function. There are lots of ways to do this. If we knew more about the function you have, we might be able to make some suggestions.

The most general purpose approach is probably the incremental delaunay construction that Damian Sheehy describes in this post on Loren's blog. Basically, you would convert your polygon into a constrained delaunayTriangulation, and then keep inserting additional points and retriangulating until you have small enough triangles.

That approach might look something like this.

Starting with my original example:

x = linspace(0,2*pi,100)';
y1 = sin(x);
y2 = 1+sin(x);
xp = [x;flipud(x)];
yp = [y1;flipud(y2)];

Create a constrained delaunayTriangulation. There are two tricky bits here. We need to pass in the edges of the polygon (so no triangles cross it), and we need to throw out any triangles which are outside the polygon:

edges = 1:numel(xp);
edges(2,:) = circshift(edges(1,:),[0 -1]);
dt = delaunayTriangulation([xp, yp],edges');
tin = isInterior(dt);
triplot(dt(tin,:),dt.Points(:,1),dt.Points(:,2))

Now we need to add more vertices. You'll want to figure out a way to do this that takes advantage of what you know about the function behind your colors. I'll just use a random number generator.

dt.Points(end+(1:numel(xn)),:) = [xn', yn'];
tin = isInterior(dt);
triplot(dt(tin,:),dt.Points(:,1),dt.Points(:,2))

Now we can compute the colors for the points which are on or in the polygon:

intris = dt(tin,:);
inpoints = unique(intris(:));
xp = dt.Points(:,1);
yp = dt.Points(:,2);
cp = nan(size(xp));
cp(inpoints) = cos(xp(inpoints)) .* cos(yp(inpoints));

And then plot them with patch (which is the low-level command that fill was calling for you):

patch('Vertices',[xp,yp], ... 'Faces',dt(tin,:), ... 'FaceVertexCData',cp, ... 'FaceColor','interp','EdgeColor','none')

So that's the basic idea, but I can't help you with the step of inserting the new points because I don't know anything about the function underlying your colors, or about your polygon.

What makes this very, very strange is that the interpolation is 
darn-near perfect if I plot my X values as y's and viceversa:

I would expect that's just random luck. Sometimes the triangulation just happens to occur in a way that looks "reasonable". You can see that in some of the examples in that blog post. For example, the one where color is a function of x^2 looks pretty reasonable, but if you look closely you'll notice that the zero is offset to the right.

[Sarah]

Thank you so much for your (detailed!) solution! I am working on implementing it now, but I wanted to put it out there that it likely (as you said) an issue with R2014a-- I borrowed a computer running R2015a last night, and it produced a map that had no (or seemingly no) vertical color variation in the fill. I saved it as a .fig file and opened it on my own computer, only to see that the streaks had reappeared. It could be random luck, but it's certainly suggestive of changes in the code of fill (or patch).

[Mike Garrity]

Those two versions use different polygon triangulators, so the way they subdivide the polygon can be completely different.

But, I would caution against considering it a "fix". It's the sort of thing that can cause you a lot of heartburn. Cases where results appear reasonable, but aren't accurate can be even more dangerous than cases where you clearly get a wrong answer. In release R2014b, we strongly considered making patch error out if the color data wasn't a linear function of the coordinates, so that users would no they were in danger, but there were complaints about it erroring in cases where it "seemed to work".