Best Way to find density of [x,y] coordinates within set window

I am using Matlab 2015. What is the best way to code a set window and have it move through some [x,y] coordinate data (from start to end), count the plots inside the window for the current iteration and save those values in a seperate array? The size of the window can remain the same as it moves throughout, and each window can be centered around the coordinate point.

I am doing this to detect clusters of coordinates. Below is an illustration which hopefully aids in explaining what I am after.

[Walter Roberson]

When you move the windows through, are the x and y coordinates of the box changing by a constant difference each time and the box is remaining the same size? The example you posted uses different sizes of box and non-constant intervals, which makes it harder to vectorize.

[James Morris]

I would like the box to remain the same size, but yes the intervals between [x,y] plots are not constant.

These coordinates represent the trajectory of a moving object over time. The clusters are created when the object is stationary for a short moment of time.

I need to be able to move through the data in set windows and detect the clusters. Detecting the clusters in the most efficient way and identifying what index they are of the dataset is my ultimate goal.

[Adam Danz]

If the boxes jumped to the next data point and were always centered at a data point, the solution would become much easier. But in your diagram it seems like you're mentally interpolating the trajectory and sliding the window along the interplation. Since your x coordinates are not unique along the trajectory, you can't use matlab's interp1() or spline() or pchip(), etc.

FYI, note a previous discussion on this topic. https://www.mathworks.com/matlabcentral/answers/475869-how-do-i-measure-density-of-data-points-within-a-window-in-plots-such-as-the-following

[James Morris]

Adam Danz, having the box centered at each data point as it jumped between them would be fine for me. What do you recommend? I'm using Matlab 2015

[Adam Danz]

I added an answer based on the info we gattered here. My answer contains an optional version that animates the window and displays the number of points within the window. It may be helpful to visualize the window.

[Walter Roberson]

If the point is detecting clusters, then I do not see why the increment should be non-constant.

[Adam Danz]

To add to Walter's comment above, if the point is to detect clusters, histcounts2() or histrogram2() might be more suitable (as I mentioned in your first question).

[Walter Roberson]

To handle the moving window, you can do histcounts2() using bins that are the width of a constant increment in x and y, after which you would do conv2() to calculate the total number of entries that are in each larger box -- like a moving sum but in 2D. This can easily be vectorized -- but it does depend upon the size of the bins being known ahead of time.

Answer 1

[Adam Danz]

Here's an alternative that uses inpolygon() to determine whether coordinates are within a defined rectangular window centered at each of your coordinates. The size of the window can be adjusted in the variable "windowSize".

'inCount' is the number of points within the window defined by windowSize, centered on each coodinate.

'onCount' is the number of points that fall exactly on the window (which will be rare).

It would be very easy to also store the coordinates of all points that fall into each window if you needd that information, too.

load('xy_coordinates.mat'); % better to use a full path

windowSize = [20,20]; % window size, [width,height]; 
inCount = nan(size(xy_coordinates,1),1); 
ouCount = inCount; 
onCount = 0; 
% Loop through each data point
for i = 1:size(xy_coordinates,1)
    % Determine square window coordinates, centered at i-th datapoint
    windCoord = xy_coordinates(i,:) + windowSize/2.*[1,1;-1,-1];  %[UpperLeft, UpperRight; LowerLeft, LowerRight] coordinates
    % Alternative: windCoord = xy_coordinates([i,i],:) + [windowSize/2;-windowSize/2]
    % Determine how many points are within window
    [in,on] = inpolygon(xy_coordinates(:,1),xy_coordinates(:,2),windCoord([2,1,1,2,2]),windCoord([4,4,3,3,4])); 
    % to plot the window : hold on; plot(windCoord([2,1,1,2,2]),windCoord([4,4,3,3,4]), 'r-')
    % Count the number of points IN and ON the window
    inCount(i) = sum(in); 
    onCount(i) = sum(on); 
end

With animation....

In this version below, you can watch the animated window move through your data and there is text in the upper left corner of the plot showing the point counts within the window.

load('xy_coordinates.mat'); % better to use a full path

figure(); 
plot(xy_coordinates(:,1), xy_coordinates(:,2), '-ob')
axis equal
hold on

windowSize = [20,20]; % window size, [width,height]; 
inCount = nan(size(xy_coordinates,1),1); 
ouCount = inCount; 
onCount = 0; 
wh = gobjects(1); 
th = text(min(xlim()),max(ylim()),'Count: 0', 'VerticalAlignment', 'top','FontSize',14);
% Loop through each data point
for i = 1:size(xy_coordinates,1)
    % Determine square window coordinates, centered at i-th datapoint
    windCoord = xy_coordinates(i,:) + windowSize/2.*[1,1;-1,-1];  %[UpperLeft, UpperRight; LowerLeft, LowerRight] coordinates
    % Alternative: windCoord = xy_coordinates([i,i],:) + [windowSize/2;-windowSize/2]
    % Determine how many points are within window
    [in,on] = inpolygon(xy_coordinates(:,1),xy_coordinates(:,2),windCoord([2,1,1,2,2]),windCoord([4,4,3,3,4])); 
    % Count the number of points IN and ON the window
    inCount(i) = sum(in); 
    onCount(i) = sum(on); 
    
    % Draw window and label number of memebers
    delete(wh)
    wh = plot(windCoord([2,1,1,2,2]),windCoord([4,4,3,3,4]), 'r-');
    th.String = sprintf('Count: %d',inCount(i)); 
    drawnow
    pause(.1)
end

[James Morris]

Hi Adam Danz. I am having issues with this line:

windCoord = xy_coordinates(i,:) + windowSize/2.*[1,1;-1,-1];

xy_coordinates(i,:) + windowSize/2 returns a size of (1 2) and [1 1 ; -1 -1] is size (2 2) thus I am getting the error:

Error using .*

Matrix dimensions must agree.

Any idea how to best rectify this?

Thanks

[Adam Danz]

That operation is OK with more recent releases of matlab but not with older releases. I updated my answer with an alternative line:

windCoord = xy_coordinates([i,i],:) + [windowSize/2;-windowSize/2]; 

[James Morris]

Amazing, thank you. Loving that animation!

[Adam Danz]

Glad that worked out! You can apply other types of animations such as building the count-plot that Image Analyst in real-time. The possibilities are limitless ;)

Answer 2

[Image Analyst]

James:

Attached is how you do it using a rectangular window. Choose the window width to be square if you want. It counts the number of other data points ANYWHERE in the data. If you want it just considering the past N elements (like recent history), you'd need to change the deltax and deltay computations to just include x and y from k-N to k ONLY instead of the whole vector like this: x(max([1, k-N]) : k). Same for y. Let me know how it goes.

clc;    % Clear the command window.
close all;  %Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

% Get original data.
s = load('xy_coordinates.mat')
x = s.xy_coordinates(:, 1);
y = s.xy_coordinates(:, 2);

% Plot original data.
subplot(2, 1, 1);
plot(x, y, 'r.-', 'LineWidth', 2, 'MarkerSize', 20);
grid on;
xlabel('x', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst for James Morris', 'NumberTitle', 'Off')

% Make both x and y axes have the same scale so it's not distorted and
% we get a true idea of what the window looks like.
axis equal;
xticks(0:50:1000);

% Now find out how many are within a sliding window centered about each data point.
windowWidthX = 30;  % Adjust as needed.
windowWidthY = 50;  % Adjust as needed.
% Plot the window over the first point so we can see its size.
hold on;
rectangle('Position', [x(1)-windowWidthX/2, y(1)-windowWidthY/2, windowWidthX, windowWidthY]);
caption = sprintf('Box size: %.1f tall by %.1f wide.  One box is shown over the very first point.', windowWidthY, windowWidthX);
title(caption, 'FontSize', fontSize);

withinWindow = zeros(1, length(x)); % Initialize for speed.
for k = 1 : length(x)
	thisX = x(k);
	thisY = y(k);
	% Compute distance from ALL other data points to this data point
	% (not just to those within a certain distance in the past along the trace).
	deltax = abs(thisX - x);
	deltay = abs(thisY - y);
	indexesInsideWindow = abs(deltax <= windowWidthX/2) & abs(deltay < windowWidthY/2);
	withinWindow(k) = sum(indexesInsideWindow);
end

% Plot point density as a function of index number.
subplot(2, 1, 2);
plot(1:length(withinWindow), withinWindow, 'b.-', 'LineWidth', 2, 'MarkerSize', 20);
grid on;
xlabel('Element number (index)', 'FontSize', fontSize);
ylabel('Count', 'FontSize', fontSize);
title('Count within window.', 'FontSize', fontSize);

[James Morris]

This seems perfect. Thank you so much Image analyst. I'll let you know how it goes

[James Morris]

Image Analyst this worked great for me, although I had issues with xticks(0:50:1000). I am using Matlab 2015, so I had to change that, but apart from all that it's perfect. Thank you

[Image Analyst]

Yes, that's a newer function and it's not crucial. You can take it out and it will just automatically decide on the tick mark values.

Since it worked perfectly for you, maybe you can "Accept this Answer" and vote for any other answers that also worked (if you tried them).

[James Morris]

I wish i could accept all these answers, they've been so helpfu! I've voted for everything you've posted. Thanks Image Analyst

Answer 3

[Image Analyst]

Can you supply your x and y in a .mat file

save('answers.mat', 'x', 'y');

Use the paper clip icon.

I'll probably just scan along and for each point, find the distance to each point, then count the number less than some window (radius) that you want to consider.

for k = 1 : length(x)
    thisX = x(k);
    thisY = y(k);
    allDistances = sqrt((thisX - x).^2 + (thisY-y).^2);
    withinWindow(k) = sum(allDistances < someValue);
end

withinWindow is a vector that will contain the count of all other points that are closer than "someValue" to the point at that index.

[Walter Roberson]

Note that this will use circular windows, not rectangular windows.

[James Morris]

Thank you, I have attached a .mat file.

How would I go about using square windows? I have been trying to define the corners using:

[x(i) - diff(x), y(i) + diff(y)] (top left corner)

[x(i) + diff(x), y(i) - diff(y)] (bottom right corner)