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

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.

[Scott MacKenzie]

It might help if you post the code that generated the image in your question.

Answer 1

[Mathieu NOE]

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 :

smoothn - File Exchange - MATLAB Central (mathworks.com)

final results

a plot of the density map (smoothed) with fitted circle overlaid

and the plot of the data selection (inside circle)

code :

load('Data.mat');

[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

[Raju Kumar]

@Mathieu NOE It was very helpful. Thanks so much. I used an 'inpolygen' method to extract the data inside the circle.

[Mathieu NOE]

My pleasure

thanks forthe info about inpolygen (haven't thought about that option)