Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Matlab: patch area between two curves which depend on the curves values

I'm trying to fill an area between two curves with respect to a function which depends on the values of the curves.

Here is the code of what I've managed to do so far

    i=50;
cc = @(xx,x,y) 1./(1+(exp(-xx)/(exp(-x)-exp(-y))));
n_vec = 2:0.1:10;
x_vec = linspace(2,10,length(n_vec));
y_vec = abs(sin(n_vec));
N=[n_vec,fliplr(n_vec)];  
X=[x_vec,fliplr(y_vec)];
figure(1)
subplot(2,1,1)
hold on
plot(n_vec,x_vec,n_vec,y_vec)
hp = patch(N,X,'b')
plot([n_vec(i) n_vec(i)],[x_vec(i),y_vec(i)],'linewidth',5)
xlabel('n'); ylabel('x')
subplot(2,1,2)

xx = linspace(y_vec(i),x_vec(i),100);
plot(xx,cc(xx,y_vec(i),x_vec(i)))
xlabel('x'); ylabel('c(x)')

This code produces the following graph

enter image description here

The color code which I've added represent the color coding that each line (along the y axis at a point on the x axis) from the area between the two curves should be.

Overall, the entire area should be filled with a gradient color which depends on the values of the curves.

I've assisted the following previous questions but could not resolve a solution

MATLAB fill area between lines

Patch circle by a color gradient

Filling between two curves, according to a colormap given by a function MATLAB

NOTE: there is no importance to the functional form of the curves, I would prefer an answer which refers to two general arrays which consist the curves.

like image 384
jarhead Avatar asked Sep 04 '19 06:09

jarhead


2 Answers

The surf plot method

  1. The same as the scatter plot method, i.e. generate a point grid.
y = [x_vec(:); y_vec(:)];
resolution = [500,500];
px = linspace(min(n_vec), max(n_vec), resolution(1));
py = linspace(min(y), max(y), resolution(2));
[px, py] = meshgrid(px, py);
  1. Generate a logical array indicating whether the points are inside the polygon, but no need to extract the points:
in = inpolygon(px, py, N, X);
  1. Generate Z. The value of Z indicates the color to use for the surface plot. Hence, it is generated using the your function cc.
pz = 1./(1+(exp(-py_)/(exp(-y_vec(i))-exp(-x_vec(i)))));
pz = repmat(pz',1,resolution(2));
  1. Set Z values for points outside the area of interest to NaN so MATLAB won't plot them.
pz(~in) = nan;
  1. Generate a bounded colourmap (delete if you want to use full colour range)
% generate colormap
c = jet(100);
[s,l] = bounds(pz,'all');
s = round(s*100);
l = round(l*100);
if s ~= 0
    c(1:s,:) = [];
end
if l ~= 100
    c(l:100,:) = [];
end
  1. Finally, plot.
figure;
colormap(jet)
surf(px,py,pz,'edgecolor','none');
view(2) % x-y view

enter image description here

Feel free to turn the image arround to see how it looks like in the Z-dimention - beautiful :)

enter image description here


Full code to test:

i=50;
cc = @(xx,x,y) 1./(1+(exp(-xx)/(exp(-x)-exp(-y))));
n_vec = 2:0.1:10;
x_vec = linspace(2,10,length(n_vec));
y_vec = abs(sin(n_vec));

% generate grid
y = [x_vec(:); y_vec(:)];
resolution = [500,500];
px_ = linspace(min(n_vec), max(n_vec), resolution(1));
py_ = linspace(min(y), max(y), resolution(2));
[px, py] = meshgrid(px_, py_);

% extract points
in = inpolygon(px, py, N, X);

% generate z
pz = 1./(1+(exp(-py_)/(exp(-y_vec(i))-exp(-x_vec(i)))));
pz = repmat(pz',1,resolution(2));
pz(~in) = nan;

% generate colormap
c = jet(100);
[s,l] = bounds(pz,'all');
s = round(s*100);
l = round(l*100);
if s ~= 0
    c(1:s,:) = [];
end
if l ~= 100
    c(l:100,:) = [];
end

% plot
figure;
colormap(c)
surf(px,py,pz,'edgecolor','none');
view(2)
like image 97
Anthony Avatar answered Oct 11 '22 14:10

Anthony


You can use imagesc and meshgrids. See comments in the code to understand what's going on.

Downsample your data

% your initial upper and lower boundaries
n_vec_long = linspace(2,10,1000000);

f_ub_vec_long = linspace(2, 10, length(n_vec_long));
f_lb_vec_long = abs(sin(n_vec_long));

% downsample
n_vec = linspace(n_vec_long(1), n_vec_long(end), 1000); % for example, only 1000 points

% get upper and lower boundary values for n_vec
f_ub_vec = interp1(n_vec_long, f_ub_vec_long, n_vec); 
f_lb_vec = interp1(n_vec_long, f_lb_vec_long, n_vec); 

% x_vec for the color function
x_vec = 0:0.01:10;

Plot the data

% create a 2D matrix with N and X position
[N, X] = meshgrid(n_vec, x_vec);

% evaluate the upper and lower boundary functions at n_vec
% can be any function at n you want (not tested for crossing boundaries though...)
f_ub_vec = linspace(2, 10, length(n_vec));
f_lb_vec = abs(sin(n_vec));

% make these row vectors into matrices, to create a boolean mask
F_UB = repmat(f_ub_vec, [size(N, 1) 1]);
F_LB = repmat(f_lb_vec, [size(N, 1) 1]);

% create a mask based on the upper and lower boundary functions
mask = true(size(N));
mask(X > F_UB | X < F_LB) = false;

% create data matrix
Z = NaN(size(N));

% create function that evaluates the color profile for each defined value 
% in the vectors with the lower and upper bounds
zc = @(X, ub, lb) 1 ./ (1 + (exp(-X) ./ (exp(-ub) - exp(-lb))));
CData = zc(X, f_lb_vec, f_ub_vec); % create the c(x) at all X

% put the CData in Z, but only between the lower and upper bound.
Z(mask) = CData(mask);

% normalize Z along 1st dim
Z = normalize(Z, 1, 'range'); % get all values between 0 and 1 for colorbar

% draw a figure!
figure(1); clf;
ax = axes;                          % create some axes
sc = imagesc(ax, n_vec, x_vec, Z);  % plot the data
ax.YDir = 'normal' % set the YDir to normal again, imagesc reverses it by default;

xlabel('n')
ylabel('x')

enter image description here

This already looks kinda like what you want, but let's get rid of the blue area outside the boundaries. This can be done by creating an 'alpha mask', i.e. set the alpha value for all pixels outside the previously defined mask to 0:

figure(2); clf;
ax = axes;                          % create some axes
hold on;
sc = imagesc(ax, n_vec, x_vec, Z);  % plot the data
ax.YDir = 'normal' % set the YDir to normal again, imagesc reverses it by default;

% set a colormap 
colormap(flip(hsv(100)))

% set alpha for points outside mask
Calpha = ones(size(N));
Calpha(~mask) = 0;
sc.AlphaData = Calpha;

% plot the other lines
plot(n_vec, f_ub_vec, 'k', n_vec, f_lb_vec, 'k' ,'linewidth', 1)

% set axis limits
xlim([min(n_vec), max(n_vec)])
ylim([min(x_vec), max(x_vec)])

![enter image description here

like image 42
rinkert Avatar answered Oct 11 '22 14:10

rinkert