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

3 Ansichten (letzte 30 Tage)
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.

Akzeptierte Antwort

Mathieu NOE
Mathieu NOE am 2 Mai 2023
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'));
% Reshape nPerBin to grid size
Z = reshape(Z, [nBins, nBins]);
%% smooth it
s_factor = 10; % smoothing parameter
Zs = smoothn(Z,s_factor); % FEX :
% 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)
xcircle = R*cos(theta')+xc;
ycircle = R*sin(theta')+yc;
% plot data
ih = imagesc(XEdge, YEdge, Zs);
hold on
axis equal;
% Reverse y axis
set(gca, 'YDir', 'normal');
% Change colormap
colormap jet
% Label the axes
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
axis equal;
% Label the axes
function [xc,yc,R] = Landau_Smith(x,y)
%------This function can be used to fit circle from arc points or circle
%------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);
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
  2 Kommentare
Raju Kumar
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
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.

Weitere Antworten (0)




Community Treasure Hunt

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

Start Hunting!

Translated by