# Get quiver plot arrows to connect up contour lines

72 visualizzazioni (ultimi 30 giorni)
David Gillcrist il 20 Feb 2024
Modificato: David Gillcrist il 21 Feb 2024
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:
1. Get the arrows to start at a contour and end at the next contour. So there should be no arrows "crossing" a contour line
2. 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
f = function_handle with value:
@(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));
% 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, ...
end
% Set properties for the main axes
axis equal
xlim(ax, [-1 1])
ylim(ax, [-1 1])
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 = 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;
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 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Risposte (2)

Cris LaPierre il 21 Feb 2024
Modificato: Cris LaPierre il 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])
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.
##### 1 CommentoMostra -1 commenti meno recentiNascondi -1 commenti meno recenti
Cris LaPierre il 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));
% 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, ...
end
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

Accedi per commentare.

Adam Danz il 21 Feb 2024
Modificato: Adam Danz il 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));
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));
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)));
end
end
cmap = spring(255);
ax.Colormap = cmap;
You can even add two colorbars by following this doc example.
##### 0 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Categorie

Scopri di piÃ¹ su Red in Help Center e File Exchange

R2022a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by