# How to separate high and low concentration areas in a given dataset as below?

1 Ansicht (letzte 30 Tage)
Raju Kumar am 27 Apr. 2023
Kommentiert: Mathieu NOE am 10 Mai 2023
I have got a matrix as attached. Its first column consists of index for m-by-n (here 256-by-256) grid pixels and second column consists of corresponding intensity/count. After plotting using 'imagesc' it looks like the figure below. It can clearly be seen that middle region is more densed that the outer region. I would like to define a boundary (let's say using a circle as depicted below). How do I do that? I want the low density area to be masked later and thus have only highly concentrated area. The data is attached. Please suggest a way to do this. Your help will be greatly appreciated.
##### 1 KommentarKeine anzeigenKeine ausblenden
Scott MacKenzie am 1 Mai 2023
It might help if you post the code that generated the image in your question.

Melden Sie sich an, um zu kommentieren.

### Akzeptierte Antwort

Mathieu NOE am 2 Mai 2023
hello
try this
maybe not the shortest route to final destination, especially if like me you don't have access to hist3 function , so I had to figure out another home made solution to plot a "density" map of your data
from there I decided that the fitted circle should be based on z data in range 8 to 9 (my own decision, you can test with other values)
also to make a 2D surface smoothing of your density data I opted for this Fex submission :
final results
a plot of the density map (smoothed) with fitted circle overlaid and the plot of the data selection (inside circle) code :
[m,n] = size(data);
[x,y] = ind2sub([256,256],data(:,1));
z = data(:,2);
%% bin the data (Hist3 clone)
nBins = 50; %number of bins (there will be nBins + 1 edges)
XEdge = linspace(min(x),max(x),nBins);
YEdge = linspace(min(y),max(y),nBins);
[~, xBin] = histc(x, XEdge);
[~, yBin] = histc(y, YEdge);
% count number of elements per (x,y) pair
[xIdx, yIdx] = meshgrid(1:nBins, 1:nBins);
xyPairs = [xIdx(:), yIdx(:)];
Z = zeros(size(xyPairs,1),1);
for i = 1:size(xyPairs,1)
Z(i) = sum(ismember([xBin, yBin], xyPairs(i,:), 'rows'));
end
% Reshape nPerBin to grid size
Z = reshape(Z, [nBins, nBins]);
%% smooth it
s_factor = 10; % smoothing parameter
Zs = smoothn(Z,s_factor); % FEX : https://fr.mathworks.com/matlabcentral/fileexchange/25634-smoothn/
% select data BELOW and ABOVE thresholds
ind = find(Zs<9 & Zs>8);
[xind,yind] = ind2sub(size(Z),ind);
Xsel = XEdge(xind);
Ysel = YEdge(yind);
Zsel = Zs(ind);
% circle fit
[xc,yc,R] = Landau_Smith(Xsel,Ysel)
theta=0:pi/180:2*pi;
xcircle = R*cos(theta')+xc;
ycircle = R*sin(theta')+yc;
% plot data
figure(1)
ih = imagesc(XEdge, YEdge, Zs);
hold on
plot(xcircle,ycircle,'LineWidth',2);
axis equal;
colorbar('vert');
% Reverse y axis
set(gca, 'YDir', 'normal');
% Change colormap
colormap jet
% Label the axes
xlabel('x')
ylabel('y')
title('hist3 simulation');
% lastly keep only data inside circle
[th,r] = cart2pol(x-xc,y-yc); % convert all data to polar
ind = (r<=R);
r = r(ind);
th = th(ind);
[x,y] = pol2cart(th,r); % convert back to cartesian and add also circle center coordinates
x = x + xc;
y = y + yc;
% plot selected data
figure(2)
plot(x,y,'.',xcircle,ycircle,'LineWidth',2);
axis equal;
% Label the axes
xlabel('x')
ylabel('y')
title('Selection');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xc,yc,R] = Landau_Smith(x,y)
%------This function can be used to fit circle from arc points or circle
%------points--
%------From your master file, call below function.
%------x and y are the coordinates of the scatter points
%------xcnew and ycnew will represent the fitted circle's center
%------Rnew is the radius, units are of the same as x and y
% [xcnew,ycnew,Rnew] = Landau_Smith(x,y);
%%%-----below code is optional(just for visualization)
% theta=0:pi/180:2*pi;
% xcircle = R*cos(theta')+xc;
% ycircle = R*sin(theta')+yc;
% plot(x,y,'.',xcircle,ycircle,'LineWidth',2);
% axis equal;
%----- Dont modify anything below this line ------
N = length(x);
p1 = 0; p2 =0; p3 =0; p4=0; p5=0; p6=0; p7=0; p8=0; p9=0;
for i=1:N
p1 = p1 + x(i);
p2 = p2 + x(i)*x(i);
p3 = p3 + x(i)*y(i);
p4 = p4 + y(i);
p5 = p5 + y(i)*y(i);
p6 = p6 + x(i)*x(i)*x(i);
p7 = p7 + x(i)*y(i)*y(i);
p8 = p8 + y(i)*y(i)*y(i);
p9 = p9 + x(i)*x(i)*y(i);
end
a1 = 2 * (p1*p1 - N*p2);
b1 = 2 * (p1*p4 - N*p3);
a2 = b1;
b2 = 2 * (p4*p4 - N*p5);
c1 = p2*p1 - N*p6 + p1*p5 - N*p7;
c2 = p2*p4 - N*p8 + p4*p5 - N*p9;
xc = (c1*b2-c2*b1)/(a1*b2-a2*b1); % returns the center along x
yc = (a1*c2-a2*c1)/(a1*b2-a2*b1); % returns the center along y
R = sqrt((p2 - 2*p1*xc + N*xc*xc + p5 - 2*p4*yc + N*yc*yc)/N); % Radius of circle
end
##### 2 Kommentare1 älteren Kommentar anzeigen1 älteren Kommentar ausblenden
Mathieu NOE am 10 Mai 2023
My pleasure

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Axis Labels finden Sie in Help Center und File Exchange

R2021a

### Community Treasure Hunt

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

Start Hunting!