pcolor plot: how to scale the smoothing differently in x and y directions

I have data of say this type (this is a simplified example):

A = [0 0 0 0 0 0
    0 0 2 0 0 0
    0 0 9 0 0 0
    0 0 1 1 0 0
    0 0 0 8 0 0
    0 0 0 1 0 0
    0 0 0 0 0 0];

If I do:

figure; pcolor(A); shading interp

then I get two distinct "peaks":

Now the thing is that I know from the property of the data (two-dimensional gas chromatography) that it is actually the same peak, it is normal that in each column the peak is moving by some pixels downward (in my simplified example here it is moving by two pixels, in the actual data it is even more). I would like to tell pcolor that distance in the x and y direction does not count the same (a sort of weighting maybe?), so that pcolor instead understands that these are one and the same "peak" and plot that as a single "peak". (This must also be directional: if the peak "goes down" in the next column (by approximately the right extent), then it should be considered the same peak. Essentially, instead of considering x and y distances, I want somehow to consider the distance at an angle)... Is there an option to do that easily?

[Jonas GROS]

excuse the poor "hand" drawing, but basically schematically this is what I want to get plotted for my data in matrix A, schematically hand-drawn.

[dpb]

A = [0 0 0   0 0 0
     0 0 2   0 0 0
     0 0 9   0 0 0
     0 0 8.5 8.5 0 0
     0 0 0   8 0 0
     0 0 0   1 0 0
     0 0 0   0 0 0];
subplot(1,2,1)
pcolor(A); shading interp

As the above shows, you'll need more resolution to build an oval object based on interpolation between the two peaks in order to draw smooth boundaries.

Simply inserting one column here

A = [0 0 0 0   0 0 0
     0 0 2 0   0 0 0
     0 0 9 0   0 0 0
     0 0 1 8.5 1 0 0
     0 0 0 0   8 0 0
     0 0 0 0   1 0 0
     0 0 0 0   0 0 0];
subplot(1,2,2)
pcolor(A); shading interp

sorta' helps but they're now disjoint without connecting rows of nonzero.

Probably attaching a real example dataset would be the next step for somebody to try something "in anger"...

[Image Analyst]

What size do you want the final, single peak to be? As large as the two -- covering both locations entirely -- or just as large as one of the single peaks (but then where would it be located)?

Answer 1

[Mathieu NOE]

hello

maybe this - I opted to overlay the two peaks on the center position (which I admit is not exactly what was asked)...

but first I needed to resample the data to have more in between pixels between the two peaks (remove if necessary for future usage)

NB : I prefered here to use peaks2 instead of islocalmax2 (I don't know why I don't like it)

available here : peaks2 - find peaks in 2D data without additional toolbox - File Exchange - MATLAB Central

result should look like :

data = [0 0 0 0 0 0
    0 0 2 0 0 0
    0 0 9 0 0 0
    0 0 1 1 0 0
    0 0 0 8 0 0
    0 0 0 1 0 0
    0 0 0 0 0 0];

[m,n] = size(data);

% resample the data with factor r 
r = 4;
xx = (1:1/r:n);
yy = (1:1/r:m);
[Xq,Yq] = meshgrid(xx,yy);
dataI = interp2(data,Xq,Yq);

%% main code below 
[m,n] = size(dataI);
%% find the two peaks location
MinPeakH = 0.5*max(dataI,[],'all');
[pks,locs_y,locs_x] = peaks2(dataI,'MinPeakHeight',MinPeakH); 

% first peak
x01 = locs_x(1);
y01 = locs_y(1);

% second peak
x02 = locs_x(2);
y02 = locs_y(2);

%% re-centering the data to the mid point
% mid point between two peaks
xm = round(mean(locs_x));
ym = round(mean(locs_y));

% shift peak1 data 
radius_pixels = r;
data1s = ones(m,n)*min(data,[],'all');
data1s(ym-radius_pixels:ym+radius_pixels,xm-radius_pixels:xm+radius_pixels) = dataI(y01-radius_pixels:y01+radius_pixels,x01-radius_pixels:x01+radius_pixels);

% shift peak2 data 
data2s = ones(m,n)*min(data,[],'all');
data2s(ym-radius_pixels:ym+radius_pixels,xm-radius_pixels:xm+radius_pixels) = dataI(y02-radius_pixels:y02+radius_pixels,x02-radius_pixels:x02+radius_pixels);

% overlay = mean of both shifted datas
datas = (data1s+data2s)/2;

figure
subplot 121
imagesc(xx,yy,dataI); hold on
plot(xx(locs_x),yy(locs_y),'+r','markersize',25);
plot(xx(xm),yy(ym),'+k','markersize',25);
subplot 122
imagesc(xx,yy,datas);  hold on
plot(xm,ym,'+k','markersize',25);

[Mathieu NOE]

not sure if this is a better approach - tried somehow to fill (interpolate) the area between the two peaks

NB that to get to this visual effect I needed quite a bit of tweaking my home made interpolation and also needed a bit of 2D smoothing to get the final touch

for the smoothing I used : smooth2a - File Exchange - MATLAB Central

you can try with your own prefered solution if you want

Nb that smoothing will cause the peaks amplitude to be reduced (or in other words, the peak width increase & amplitude decreases at the same time)

without smoothing the result is not as nice looking .

my result so far :

code :

data = [0 0 0 0 0 0
    0 0 2 0 0 0
    0 0 9 0 0 0
    0 0 1 1 0 0
    0 0 0 8 0 0
    0 0 0 1 0 0
    0 0 0 0 0 0];

[m,n] = size(data);

% resample the data with factor r 
r = 4;
xx = (1:1/r:n);
yy = (1:1/r:m);
[Xq,Yq] = meshgrid(xx,yy);
dataI = interp2(data,Xq,Yq);

%% main code below 
[m,n] = size(dataI);
%% find the two peaks location
MinPeakH = 0.5*max(dataI,[],'all');
[pks,locs_y,locs_x] = peaks2(dataI,'MinPeakHeight',MinPeakH); 

% first peak
x01 = locs_x(1);
y01 = locs_y(1);

% second peak
x02 = locs_x(2);
y02 = locs_y(2);

%fill the area between the peaks (along vertical direction)
xx = min(locs_x):max(locs_x);
yy = min(locs_y):max(locs_y);

ytop = dataI(min(locs_y),xx);
ybottom = dataI(max(locs_y),xx);

data_out = dataI;

for k = 1:numel(yy)
    a = (k-1)/(numel(yy)-1);
    ymix = (1-a)*ytop + a*ybottom ;
    % amplify peak amplitude (optionnal)
    avg_peak = (1-a)*max(ytop) + a*max(ybottom) ; % linear interpolation
    [val,~] = max(ymix);
    corr_factor = avg_peak/val;
    corr_factor = 1 + 0.25*(corr_factor-1);
    data_out(yy(k),xx) = corr_factor*ymix;

end

% smooth it 
data_out = smooth2a(data_out,2,2);

figure
subplot(1,2,1)
imagesc(dataI); hold on
plot(locs_x,locs_y,'+r','markersize',25);
colorbar('vert')
subplot(1,2,2)
imagesc(data_out); hold on
plot(locs_x,locs_y,'+r','markersize',25);
colorbar('vert')

[Mathieu NOE]

yet another try - somehow a bit inspired by @Star Strider answer using contour

here I wanted to avoid the narrowing in the center area so using boundary and fill to get this result :

A = [0 0 0 0 0 0
    0 0 2 0 0 0
    0 0 9 0 0 0
    0 0 1 1 0 0
    0 0 0 8 0 0
    0 0 0 1 0 0
    0 0 0 0 0 0];

%% main code  
figure
subplot(1,2,1)
imagesc(A);
set(gca,'YDir','normal')
colorbar('vert')
subplot(1,2,2)
hold on
cmap = colormap;
nmap = size(cmap,1);

for level = min(A(:)):max(A(:))

    [C,h] = contour(A,level*[1 1],'edgecolor','none');
    ind = find(C(1,:)~=level); 
    % get all x,y coordinates for this level
    xc = C(1,ind);
    yc = C(2,ind);
    k = boundary(xc(:),yc(:),0);
    % fill the area (with the right color)
    ind_map = round((level - min(A(:)))/(max(A(:))-min(A(:))) * (nmap-1) + 1);
    fill(xc(k),yc(k),cmap(ind_map,:),'edgecolor','none')
end
colorbar('vert')

Answer 2

[Star Strider]

I am not certain that what you want to do is possible.

The closest I can get is slightly to interpolate your data and then use the contourf function.

A = [0 0 0 0 0 0
    0 0 2 0 0 0
    0 0 9 0 0 0
    0 0 1 1 0 0
    0 0 0 8 0 0
    0 0 0 1 0 0
    0 0 0 0 0 0];

x = 1:size(A,1);
y = 1:size(A,2);
[X,Y] = ndgrid(x,y);

figure
surf(X,Y,A)
view(0,90)
title("Original 'surf' Plot")

Afcn = scatteredInterpolant(X(:), Y(:), A(:), 'natural');

sf = 2;
xr = linspace(min(x), max(x), numel(x)*sf);
yr = linspace(min(y), max(y), numel(y)*sf);
[Xr,Yr] = ndgrid(xr, yr);
Ar = Afcn(Xr, Yr);

figure
surf(Xr, Yr, Ar, EdgeColor='none')
title('Interpolated ''surf'' Plot')

% view(0,90)

figure
contourf(Xr, Yr, Ar, EdgeColor='none')
title('Interpolated ''contourf'' Plot')

Other options would require changing the data, and that is generally not appropriate.

.

Answer 3

[dpb]

"... instead of considering x and y distances, I want somehow to consider the distance at an angle"

A = [0 0 0 0 0 0
    0 0 2 0 0 0
    0 0 9 0 0 0
    0 0 1 1 0 0
    0 0 0 8 0 0
    0 0 0 1 0 0
    0 0 0 0 0 0];

Npk=2;                              % number of peaks to find
[pk,ipk]=maxk(A(:),Npk);            % locate them
[r,c]=ind2sub(size(A),ipk);         % get row, column
[r c]

ans = 2×2

3 3 5 4

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

dr=diff(flip(r)); dc=-diff(flip(c));    % compute distance from upper to lower
[dr dc]

ans = 1×2

-2 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

d=hypot(dr,dc)                          % how far did it move?

d = 2.2361

[th,rho]=cart2pol(dr,dc)                % what is the distance, direction

th = 2.6779

rho = 2.2361

The above quickie is dependent upon the first peak always being to the left or above the second; if it were possible for the second to move to the left since the search is by column, it would be found first and the distance and angle would have to make sure the one being subtracted was the lower instead of just reversing the direction as above.

I still think it imperative to have a real data set in order to be able to do anything truly useful.

Answer 4

[William Rose]

@Jonas GROS,

I recommend you filter the image with a 2-D smoothing function, or kernel. The user specifies the horizontal and vertical separation of points that should be merged. The code does the rest.

First, make a sample image with two pairs of dots.

% Make sample image
im1=zeros([90,192]);
im1(37:40,60:63)=ones(4); im1(40:43,133:136)=ones(4);
im1(49:52,66:69)=ones(4); im1(46:49,121:124)=ones(4);
figure
subplot(211), imshow(im1); title('Original')  % display image

Specify the vertical and horizontal separation of points that should be merged, in pixels. If you want to merge a pair that is oriented upper left to lower right, make h and v positive. If you want to merge a pair that is oriented lower left to upper right, make h or v negative. In this case, we want to merge the pair of dots on the left, which have separation v=12, h=6.

% Specify the vertical and horizontal separation of pixels that should be merged.
v=12; h=6;  % distance in vertical and horizontal pixels

Compute the kernel, which will be used for filtering. The kernel long axis is three times the short axis (3:1). If you want the kernel relative dimensions to be 2:1, or 4:1, then try

Z=max(1-(Xr/(d/1.5)).^4-(Yr/(d/3)).^4,0);

or

Z=max(1-(Xr/(d/1.5)).^4-(Yr/(d/6)).^4,0);

instead of the corresponding line below. And so on.

% Compute the kernel
theta=atan2(h,v);
d=round(sqrt(v^2+h^2));
x=-d:d; y=-d:d;
[X,Y]=meshgrid(x,y);
Xr=X*cos(theta)+Y*sin(theta);
Yr=-X*sin(theta)+Y*cos(theta);
Z=max(1-(Xr/(d/1.5)).^4-(Yr/(d/4.5)).^4,0);

Filter the image using the kernel, and plot the filtered image.

% Use kernel to filter image
im2=imfilter(im1,Z');
% Plot filtered image
subplot(212), imshow(im2); title('Filtered')

Plot the kernel function. H and V correspond to pixel values.

% Plot kernel
figure; surf(X,Y,Z)
axis equal; view(90,90)
title(['Smoothing Filter: V=',num2str(v),', H=',num2str(h)])
colorbar; xlabel('V=vert'); ylabel('H=horiz')

●

Good luck with your 2D gas chromatography analysis.