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

4 Ansichten (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 Kommentar-1 ältere Kommentare anzeigen-1 ältere Kommentare 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 KommentareKeine anzeigenKeine ausblenden
Raju Kumar am 10 Mai 2023
@Mathieu NOE It was very helpful. Thanks so much. I used an 'inpolygen' method to extract the data inside the circle.
Mathieu NOE am 10 Mai 2023
My pleasure
thanks forthe info about inpolygen (haven't thought about that option)

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Data Distribution Plots 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!

Translated by