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

2 次查看（过去 30 天）
Raju Kumar
Raju Kumar 2023-4-27
评论： Mathieu NOE 2023-5-10
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.

请先登录，再回答此问题。

采纳的回答

Mathieu NOE
Mathieu NOE 2023-5-2
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Surface and Mesh Plots 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by