Get quiver plot arrows to connect up contour lines
42 views (last 30 days)
Show older comments
I have the following code here that plots a contour plot for a given function and also plots a color-indicated directional field of the gradient on top of it. The arrows are all the same length but the difference of color indicates the magnitude. I am trying to do two things here:
- Get the arrows to start at a contour and end at the next contour. So there should be no arrows "crossing" a contour line
- Get the colormap of the contours to be different than the colormap of the directional field
I think in order to do (1) I need to get the coordinates of the contours themselves so to provide the starting points of the arrows, but I don't know how to do that. As for (2), I better there is clever way of doing this, but my only thought is to create two entirely different plots and draw them on top of each other which seems too brute force. Any help is appreciated, code is below:
clf
f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1
[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);
Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));
% Create axes
ax = axes('Color', [1 1 1]);
hold(ax,'on')
contour(ax, X,Y,Z,10)
% Demo data
scale_size = 0.25;
[X,Y] = meshgrid(-1:0.2:1,-1:0.2:1);
Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));
[U,V] = gradient(Z,0.2,0.2);
% Define the colormap for the quiver arrows.
% cmap can have any number of rows.
cmap = spring(255);
ax.Colormap = cmap;
% Assign colors based on magnitude of vectors
vectorMagnitude = hypot(U(:),V(:));
% Scale magnitudes to rows of colormap
vecMagNorm = (vectorMagnitude-min(vectorMagnitude))./range(vectorMagnitude);
vecColorIdx = round(vecMagNorm * (size(cmap,1)-1)) + 1;
% Plot the quiver data
dU = U./sqrt(U.^2 + V.^2);
dV = V./sqrt(U.^2 + V.^2);
for i = 1:numel(Z)
quiver(ax, X(i),Y(i),dU(i),dV(i), 0.1, ...
'Color', cmap(vecColorIdx(i),:), 'LineWidth',0.2,'MaxHeadSize',2)
end
% Set properties for the main axes
axis equal
xlim(ax, [-1 1])
ylim(ax, [-1 1])
% Add colorbar
cb = colorbar(ax);
% Set colorbar range
caxis(ax, [floor(vectorMagnitude(1)), ceil(vectorMagnitude(2))])
% Label the colorbars
ylabel(cb,'Vector magnitude')
hold off
EDIT:
I was able to produce something for the vectors being spaced exactly between contours pertaining to (1) which uses the File Exchange function intersections found here: Fast and Robust Curve Intersections. The code and plot are below
clf
f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1;
[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);
Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));
syms x y
dz = gradient(f(x*cos(pi/6)-y*sin(pi/6),y*cos(pi/6)+x*sin(pi/6)),[x,y]);
dz = matlabFunction(dz);
% Create axes
ax = axes('Color', [1 1 1]);
hold(ax,'on')
NumOfContours = 10;
[M,~] = contour(ax, X,Y,Z,NumOfContours);
i = 1;
Ns = zeros(1,NumOfContours);
counter = 0;
while i<=length(M)
counter = counter + 1;
Ns(counter) = M(2,i);
i = i + M(2,i)+1;
end
X = nan(1,1000); % Pre-allocation
Y = X;
Ms = zeros(1,NumOfContours-1);
startIdx = 2;
overallCounter = 1;
for j=1:NumOfContours-1
subX = M(1,startIdx:(startIdx + Ns(j) - 1));
subY = M(2,startIdx:(startIdx + Ns(j) - 1));
startIdx = startIdx + Ns(j) + 1;
dX = gradient(subX);
dY = gradient(subY);
if j==1
PrimaryLength = sum(hypot(dX,dY));
NumVecPerUnit = 4;
dl = PrimaryLength/(NumVecPerUnit+1);
end
Lens = cumsum(hypot(dX,dY));
counter = 1;
gatePass = 1;
while gatePass
[~,idx] = min(abs(dl*counter - Lens));
if idx+1 > length(Lens)
gatePass = 0;
Ms(j) = counter-1;
else
counter = counter + 1;
X(overallCounter) = subX(idx);
Y(overallCounter) = subY(idx);
overallCounter = overallCounter + 1;
end
end
end
X(isnan(X)) = [];
Y(isnan(Y)) = [];
U = zeros(1,length(X));
V = U;
for i=1:length(X)
d = dz(X(i),Y(i));
U(i) = d(1);
V(i) = d(2);
end
% Define the colormap for the quiver arrows.
% cmap can have any number of rows.
cmap = parula(255);
ax.Colormap = cmap;
scale_size = zeros(1,length(X)); % Pre-allocate scaling array
dU = U./sqrt(U.^2 + V.^2);
dV = V./sqrt(U.^2 + V.^2);
startIdx = Ns(1)+3;
j = 0;
for i=1:length(X)
j_gate = i > cumsum([0,Ms]);
j_switch = find(j_gate == max(j_gate),1,"last");
if j_switch > j
j = j_switch;
x1 = M(1,startIdx:(startIdx + Ns(j+1) - 1));
y1 = M(2,startIdx:(startIdx + Ns(j+1) - 1));
startIdx = startIdx + Ns(j+1) + 1;
end
x2 = [X(i),X(i)+sqrt(2)*dU(i)];
y2 = [Y(i),Y(i)+sqrt(2)*dV(i)];
[x0,y0] = intersections(x1,y1,x2,y2,1);
false_switch = (isempty(x0) || isempty(y0));
x0(false_switch) = inf;
y0(false_switch) = inf;
scale_size(i) = hypot(X(i)-x0(1),Y(i)-y0(1));
end
dU = scale_size.*dU;
dV = scale_size.*dV;
remove_id = false(1,length(X));
for i=1:length(X)
if norm([X(i)+dU(i),Y(i)+dV(i)],"inf") > 1
remove_id(i) = true;
end
end
X(remove_id) = [];
Y(remove_id) = [];
dU(remove_id) = [];
dV(remove_id) = [];
quiver(ax, X+0.05*dU,Y+0.05*dV,0.925*dU,0.925*dV, ...
'Color', [0.1 0.1 0.1], 'LineWidth',0.8,'MaxHeadSize',2,'AutoScale','off')
% scatter(X,Y)
% Set properties for the main axes
axis equal
xlim(ax, [-1 1])
ylim(ax, [-1 1])
hold off
0 Comments
Answers (2)
Cris LaPierre
on 21 Feb 2024
Edited: Cris LaPierre
on 21 Feb 2024
Getting the coordinates of your contour lines is fairly straightforward, as it's an output of the contour function. There's a little coding work needed to turn that output into something you can use, but doable.
Getting the quiver arrows to appear as you describe is not giong to be an easy task. Doable, but not pretty, as that is not what a quiver plot does. I by no means attempted to arriave at a final solution, but here is a working example (simplified)
f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1;
[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);
Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));
% Create axes
ax = axes('Color', [1 1 1]);
hold(ax,'on')
%% Extract contour line coordinates
m = contour(ax, X,Y,Z,10);
i = 1;
lvl = [];
x = [];
y = [];
while i < length(m)
clvl = m(1,i);
num = m(2,i);
lvl = [lvl;clvl];
% extracted (x,y) coordinates
tmpx = m(1,i+1:i+num);
tmpy = m(2,i+1:i+num);
% convert contour coordinates to 10 equally spaced coordinates
newx = linspace(tmpx(1),tmpx(end),10);
newy = interp1(tmpx,tmpy,newx);
x = [x,newx'];
y = [y,newy'];
i = i+num+1;
end
% Added this to visualize the points
scatter(x,y,10,"magenta","filled")
%% Calculate U as difference of x coordinates and V as differeing of y coordinates
U = diff(x,[],2);
V = diff(y,[],2);
% Create quiver plot
quiver(ax,x(:,1:end-1),y(:,1:end-1),U,V,"off")
% Set properties for the main axes
axis equal
xlim(ax, [-1 1])
ylim(ax, [-1 1])
% Add colorbar
cb = colorbar(ax);
hold off
As for colormap, that is an axiss property. The solution is to either manually set the color of one of your plots, or overlay 2 axes, and specify the colormap of each. Place your plots on the axes whose colormap you want it to use.
Here are some posts about this
1 Comment
Cris LaPierre
on 21 Feb 2024
Here's a different example that shows how you could use 2 different colormaps in a single figure. This is meant to address question 2, and does not make an attempt at incorporating a solution for your first question.
f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1;
[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);
Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));
% Create axes
t = tiledlayout(1,1);
ax = axes(t,'Color', [1 1 1]);
contour(ax,X,Y,Z,10)
% Set properties for the main axes
axis equal
xlim(ax, [-1 1])
ylim(ax, [-1 1])
% Demo data
scale_size = 0.25;
[X,Y] = meshgrid(-1:0.2:1,-1:0.2:1);
Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));
[U,V] = gradient(Z,0.2,0.2);
% Define the colormap for the quiver arrows.
% cmap can have any number of rows.
ax2 = axes(t);
colormap(ax2,'spring');
cmap = colormap(ax2);
hold on
% Assign colors based on magnitude of vectors
vectorMagnitude = hypot(U(:),V(:));
% Scale magnitudes to rows of colormap
vecMagNorm = (vectorMagnitude-min(vectorMagnitude))./range(vectorMagnitude);
vecColorIdx = round(vecMagNorm * (size(cmap,1)-1)) + 1;
% Plot the quiver data
dU = U./sqrt(U.^2 + V.^2);
dV = V./sqrt(U.^2 + V.^2);
for i = 1:numel(Z)
quiver(ax2, X(i),Y(i),dU(i),dV(i), 0.1, ...
'Color', cmap(vecColorIdx(i),:), 'LineWidth',0.2,'MaxHeadSize',2)
end
% Add colorbar
cb = colorbar(ax2);
% Set colorbar range
caxis(ax2, [floor(vectorMagnitude(1)), ceil(vectorMagnitude(2))])
% Label the colorbars
ylabel(cb,'Vector magnitude')
% Set properties for the main axes
axis equal
xlim(ax2, [-1 1])
ylim(ax2, [-1 1])
ax2.Color = 'none';
hold off
Adam Danz
on 21 Feb 2024
Edited: Adam Danz
on 21 Feb 2024
Let's break this down into steps. This solution doesn't address all of the requirements in the question.
Get the arrows to start at a contour
Step 1 is getting the contour line coordinates. This isn't easy so there are several functions on the file exchange that solve this problem. I'll use getContourLineCoordinates which you'll need to download and add to the MATLAB path.
The contour lines coordinates are not sampled at a uniform interval but you probably don't want an arrow at each coordinate anyway. You could resample each curve to get coordinates along the curves at a uniform interval but I used an easier approach. The variable qint=3 means use every 3 points along the curve.
To get the horizontal and vertical gradient components at each of those coordinates, use interp2 to interpolate the gradient matrices. Then use the result to plot the quiver arrows.
I also recommending turning auto scaling off (5th input in quiver, below).
f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1;
[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);
Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));
[U,V] = gradient(Z);
ax = axes();
hold(ax,'on')
[~, h] = contour(ax, X,Y,Z,10);
% https://www.mathworks.com/matlabcentral/fileexchange/74010-getcontourlinecoordinates
CT = getContourLineCoordinates(h);
qint = 3; % quiver interval along contour lines (non-zero integer)
for i = unique(CT.Group)'
cx = CT.X(CT.Group == i);
cy = CT.Y(CT.Group == i);
idx = mod(1:numel(cx),qint)==0;
uInterp = interp2(X,Y,U,cx(idx),cy(idx));
vInterp = interp2(X,Y,V,cx(idx),cy(idx));
quiver(cx(idx),cy(idx),uInterp,vInterp,'off','k')
end
cmap = spring(255);
ax.Colormap = cmap;
Get the arrows to start at a contour and end at the next contour
🤔 this is a head scratcher. Closer spacing of contours indicate that elevation is changing faster over a smaller distance. In other words, there's a steeper gradient. A steeper gradient would be depicted as longer quiver arrows. So, there's an inverse releationship between the distance between contours and the length of quiver arrows. I'll skip this step.
So there should be no arrows "crossing" a contour line
If this is what you're concerned with, you could try scaling the quiver magnitudes using the 5th input argument that I currently have set to 'off', although it may be difficult to algorithmically choose a scale factor.
Color the quiver arrows based on magnitude
This we can do. But it will require plotting each arrow individually. In this case, it's critical to turn off auto-scaling so we can be sure maintain the relative differences between each arrow. Plotting single quiver arrows results in different arrow head sizes so this solution also increases the maxHeadSize to compensate for smaller arrow heads.
f = @(x,y) 0.5*sin(pi/2*(2*x-y)) + 2*y + 1
[X,Y] = meshgrid(-1:0.05:1,-1:0.05:1);
Z = f(X*cos(pi/6)-Y*sin(pi/6),Y*cos(pi/6)+X*sin(pi/6));
[U,V] = gradient(Z);
ax = axes();
hold(ax,'on')
[~, h] = contour(ax, X,Y,Z,10);
% https://www.mathworks.com/matlabcentral/fileexchange/74010-getcontourlinecoordinates
CT = getContourLineCoordinates(h);
UVMagnitudes = hypot(U,V);
minMag = min(UVMagnitudes,[],'all');
maxMag = max(UVMagnitudes,[],'all');
nColors = 256;
quivColors = cool(nColors);
magnitudes = linspace(minMag, maxMag, nColors)';
qint = 3; % quiver interval along contour lines (non-zero integer)
for i = unique(CT.Group)'
cxg = CT.X(CT.Group == i);
cyg = CT.Y(CT.Group == i);
idx = mod(1:numel(cxg),qint)==0;
cx = cxg(idx);
cy = cyg(idx);
uInterp = interp2(X,Y,U,cx,cy);
vInterp = interp2(X,Y,V,cx,cy);
for j = 1:sum(idx)
colr = interp1(magnitudes, quivColors, hypot(uInterp(j),vInterp(j)));
quiver(cx(j),cy(j),uInterp(j),vInterp(j),'off','k','color',colr,'MaxHeadSize',.8)
end
end
cmap = spring(255);
ax.Colormap = cmap;
0 Comments
See Also
Categories
Find more on Graphics Object Programming in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!