# Sorting and selecting points into sub matrices based on proximity to each other.

5 views (last 30 days)
Vance Blake on 19 Jun 2020
Commented: Adam Danz on 19 Jun 2020
Hi, I am sorting groups of points based on their proximity to the six points plotted on the boundary of the circle. Once the are sorted into (earlier half of code below), I am trying to select the points in a each subgroup matrix that have only one nearest neighbor a distance of 8 units away (latter half of code). I calcuate the distances between the points in the subgroup and then select the ones that have only one point in that subgroup that is 8 units away. I am having a problem with selecting just the points that have one nearestt neighbor (1, 12, 17, 18,19). I would greatly appreicate any help with this issue. Thank you for your time.
figure
();
% Axis Labels
axis([-100, 100 -100 100]);
fontsize = 12;
xlabel('X', 'Fontsize', fontsize);
ylabel('Y', 'Fontsize', fontsize);
hold on
% Plot Outerboundary Circular Plane
x_c = 0;
y_c = 0;
center_plane = [x_c, y_c]; % center point of circular plane
SizeCenter_Plane = size(center_plane,1);
CenterPlane_x = center_plane(:,1);
CenterPlane_y = center_plane(:,2);
hold on
%%
VNA = [0,80;69.2820323027551,40;69.2820323027551,-40;0,-80;-69.2820323027551,-40;-69.2820323027551,40;-0.0991889411901791,72.0006149265118;62.0092635350352,36.6673082273818;-63.2899442743164,34.6995395426962;1.74788106683108,-72.1932777828199;-62.4959472931593,-35.7633680543928;62.4229420020957,-35.8826124487242;-0.0998391324817250,64.0006149529336;54.7364947673152,33.3346164547635;-57.2978562458777,29.3990790853924;3.49576213366216,-64.3865555656398;-55.7098622835635,-31.5267361087857;55.8427494978502,-31.3327699671573;-57.8233835908142,-42.2569917305483;-0.306917232338454,56.0032954857271;48.5101288958709,28.3114333966904;-51.3057682174389,24.0986186280885;5.24364320049324,-56.5798333484597;-49.9750731423205,-25.9489055712486;49.2625569936047,-26.7829274855904;-8.09983706549959,63.9948641521246;-62.5543739803260,23.3684057144138;11.3621197247966,-65.8427211635033;59.2926597806400,-24.1148672911231;-53.1508198884690,-48.7506154067037;1.37369501047041,48.1818159268891;41.8541634034021,23.8731766034138;-43.5766929317874,26.1629333397375;2.85168633501620,-48.9457949444629;-43.9239872279001,-20.7158988584396;43.3919012899549,-21.3482792348565;-7.36645737685568,-76.8798548567465;50.5243968185633,-37.3089832145842;12.7383455385123,-53.7817202584485;-16.0998349985174,63.9891133513155;64.3938426862084,-17.9522434210547;-45.8443047451164,-52.0086572204255;4.42916682376490,40.7882990371970;35.4144641270267,19.1265701072052;-36.4649486316018,22.4991793893776;0.0696697018587238,-41.4451027462165;-39.7219315335501,-13.9083495239449;37.5212455863052,-15.9136309841226;-43.4913854881796,34.1624784918111;-24.0596597631192,63.1883701235374;-38.5377896017639,-55.2666990341473;7.22891264953506,33.2942064763139;28.9747648506513,14.3799636109965;-35.5198758392000,-7.10080018945010;31.6505898826555,-10.4789827333888;20.4163945563335,-56.0283985353575;11.4130897549582,44.6901986375394;-46.5301493027589,-9.70737692143613;41.9312472743011,-9.23891113143420;2.73824698106806,26.6734873421071;22.5350655742760,9.63335711478786;-27.5207251612244,-7.21736961529872;24.5219452703326,-14.1097438094168;64.3648003337138,29.0219536646377;-31.4817515692312,60.2029801621344;-31.0374084890469,-52.4838438162102;28.0944435741547,-58.2750768122666;18.3970126861515,48.5920982378818;-19.5215744832488,-7.33393904114735;-60.8414503911938,15.5539388242290;46.9498522520006,17.7060090883978;-29.3266608393915,18.8874142536604;-28.1094847324876,52.9484769826052;-66.0324580639326,-28.5874993565832;16.7055570530325,-46.8346897562081;-14.9189668598841,-74.2417493795593;42.8552652265820,-39.5859154983543;26.3507145241858,49.4515323957145;-12.0309445868776,-4.52494211064166;-62.0630131973175,7.64775236610042;-36.6318973496899,-24.0061002516836;-4.53735286063732,-34.9048165820909;67.2981552125186,21.5791430963785;-23.0565300842135,-53.0366366230219;-23.8025442720866,24.6739621149354;-59.1603900887728,41.5513117388991;62.5902294640123,-10.1582085645302;16.0033601290985,5.01417710341792;-67.6449233765971,-20.7516869790786;7.83895136177013,72.9935542653020;-0.826382771885093,-27.8175956175389;-21.8305993620089,47.9910987940474;-17.6628453058907,19.5452101610490;70.8169082870940,14.3945500807536;52.7102064758139,41.0737479150583;21.1007545627210,-40.1502123008772;30.5983588896414,56.2307297925022;-68.3295011127666,2.67471212842579];
VNAx = VNA(:,1);
VNAy = VNA(:,2);
SizeXY = size(VNA,1);
plot(VNAx, VNAy, '*', 'color', 'k', 'markersize', 12);
for i = 1:SizeXY
centers_nodeXY = [VNAx(i), VNAy(i)];
end
%%
% Find Branch Roots on the Edge of the Boundary
ini_branch_dist = nan(numel(VNAx));
CenterEdge_threshold = 80;
for i = 1:SizeXY
for j = 1:SizeCenter_Plane
ini_branch_dist(i,j) = sqrt((CenterPlane_x(j)-VNAx(i)).^2 + (CenterPlane_y(j)-VNAy(i)).^2);
end
end
% find the points that are equal to CenterEdge_threshold:
i2keepEdge_Center = find(ini_branch_dist == CenterEdge_threshold);
% put those into one pair of arrays
keep_x = VNAx(i2keepEdge_Center);
keep_y = VNAy(i2keepEdge_Center);
% and the others into another pair of arrays
x_remaining = VNAx;
y_remaining = VNAy;
x_remaining(i2keepEdge_Center) = [];
y_remaining(i2keepEdge_Center) = [];
BranchRoots = [keep_x, keep_y];
RemainingVNs = [x_remaining, y_remaining];
%%
% Computes distances between BranchRoots along circle edge and all remaining vein nodes
distRoot_Branch = pdist2(BranchRoots, RemainingVNs);
%%
% Split RemainingVNs into groups that are nearest to Starting Root VNs in BranchRoots
[~, minRowIdx] = min(distRoot_Branch,[],1);
[GroupID, GroupList] = findgroups(minRowIdx);
Root_BranchNeighborGroups = splitapply(@(x){x},RemainingVNs,GroupID(:));
SizeRoot_BranchNG = size(Root_BranchNeighborGroups,1);
%%
%Create Branch Matrices
Branch1 = [BranchRoots(1,:);Root_BranchNeighborGroups{1}];
Branch1x = Branch1(:,1);
Branch1y = Branch1(:,2);
SizeofBranch1 = size(Branch1,1);
%%
% Find Branch 1 Terminating VNs or Branch Ends
Branch1_dist = nan(numel(Branch1x)); % places nans on the diagonal after distance calculation
VN_VN_proximity = 8;
for i = 1:SizeofBranch1 % looping from largest index to avoid calculating the size of the elim_dist2 matrix without pying the price of dynamic growing of a matrix
for j = 1:(i-1) % distance symmetric (l(i1,i2) == l(i2,i1)) so calculate it once
Branch1_dist(i,j) = sqrt((Branch1x(i)-Branch1x(j)).^2 + (Branch1y(i)-Branch1y(j)).^2);
Branch1_dist(j,i) = Branch1_dist(i,j);
end
end
% find the points that have exactly 1 nearest neighbor
i2keepB1_B1 = find(min(Branch1_dist) >= VN_VN_proximity); % This is the line of code I am trying to figure out.
% put those into one pair of arrays
keep_x1 = Branch1x(i2keepB1_B1);
keep_y1 = Branch1y(i2keepB1_B1);
% and the others into another pair of arrays
x_close_neighbors = Branch1x;
y_close_neighbors = Branch1y;
x_close_neighbors(i2keepB1_B1) = [];
y_close_neighbors(i2keepB1_B1) = [];
Branch1_Ends = [keep_x1, keep_y1];
%%
Adam Danz on 19 Jun 2020
I'm looking at this now.
Question #1.
I'm looking at these two nested loops (copied from your question)
% Find Branch Roots on the Edge of the Boundary
ini_branch_dist = nan(numel(VNAx));
CenterEdge_threshold = 80;
for i = 1:SizeXY
for j = 1:SizeCenter_Plane
ini_branch_dist(i,j) = sqrt((CenterPlane_x(j)-VNAx(i)).^2 + (CenterPlane_y(j)-VNAy(i)).^2);
end
end
CenterPlane_x is just a single value and SizeCenter_Plane == 1. Is that supposed to be the case? If so, why is ini_branch_dist initialized as a 98x98 matrix? Only the first column will contain non-nan values. If you step through that section of the code you'll see what I'm talking about.
Maybe this
ini_branch_dist = nan(numel(VNAx));
should be one of these instead?
ini_branch_dist = nan(1, numel(VNAx));
% --OR--
ini_branch_dist = nan(numel(VNAx), 1);

Adam Danz on 19 Jun 2020
Edited: Adam Danz on 19 Jun 2020
The 3 lines of code at the bottom show the solution. I've included all of the exporatory plots and other lines of code that I added just so you can see how other people may approach the problem and get to know that data. Anything I've added or changed has been marked with a comment including a hashtag and my initials, so search for them to see what I've done ( #AD ).
The endpoints with only 1 nearnest neightbor at approximately 8 units of distance are framed in red squares in the image below.
figure();
% Axis Labels
axis([-100, 100 -100 100]);
ax = gca(); % #AD: store handle
fontsize = 12;
xlabel('X', 'Fontsize', fontsize);
ylabel('Y', 'Fontsize', fontsize);
hold on
% Plot Outerboundary Circular Plane
x_c = 0;
y_c = 0;
center_plane = [x_c, y_c]; % center point of circular plane
SizeCenter_Plane = size(center_plane,1);
CenterPlane_x = center_plane(:,1);
CenterPlane_y = center_plane(:,2);
axis equal % #AD: make circle appear as circle
hold on
%%
VNA = [0,80;69.2820323027551,40;69.2820323027551,-40;0,-80;-69.2820323027551,-40;-69.2820323027551,40;-0.0991889411901791,72.0006149265118;62.0092635350352,36.6673082273818;-63.2899442743164,34.6995395426962;1.74788106683108,-72.1932777828199;-62.4959472931593,-35.7633680543928;62.4229420020957,-35.8826124487242;-0.0998391324817250,64.0006149529336;54.7364947673152,33.3346164547635;-57.2978562458777,29.3990790853924;3.49576213366216,-64.3865555656398;-55.7098622835635,-31.5267361087857;55.8427494978502,-31.3327699671573;-57.8233835908142,-42.2569917305483;-0.306917232338454,56.0032954857271;48.5101288958709,28.3114333966904;-51.3057682174389,24.0986186280885;5.24364320049324,-56.5798333484597;-49.9750731423205,-25.9489055712486;49.2625569936047,-26.7829274855904;-8.09983706549959,63.9948641521246;-62.5543739803260,23.3684057144138;11.3621197247966,-65.8427211635033;59.2926597806400,-24.1148672911231;-53.1508198884690,-48.7506154067037;1.37369501047041,48.1818159268891;41.8541634034021,23.8731766034138;-43.5766929317874,26.1629333397375;2.85168633501620,-48.9457949444629;-43.9239872279001,-20.7158988584396;43.3919012899549,-21.3482792348565;-7.36645737685568,-76.8798548567465;50.5243968185633,-37.3089832145842;12.7383455385123,-53.7817202584485;-16.0998349985174,63.9891133513155;64.3938426862084,-17.9522434210547;-45.8443047451164,-52.0086572204255;4.42916682376490,40.7882990371970;35.4144641270267,19.1265701072052;-36.4649486316018,22.4991793893776;0.0696697018587238,-41.4451027462165;-39.7219315335501,-13.9083495239449;37.5212455863052,-15.9136309841226;-43.4913854881796,34.1624784918111;-24.0596597631192,63.1883701235374;-38.5377896017639,-55.2666990341473;7.22891264953506,33.2942064763139;28.9747648506513,14.3799636109965;-35.5198758392000,-7.10080018945010;31.6505898826555,-10.4789827333888;20.4163945563335,-56.0283985353575;11.4130897549582,44.6901986375394;-46.5301493027589,-9.70737692143613;41.9312472743011,-9.23891113143420;2.73824698106806,26.6734873421071;22.5350655742760,9.63335711478786;-27.5207251612244,-7.21736961529872;24.5219452703326,-14.1097438094168;64.3648003337138,29.0219536646377;-31.4817515692312,60.2029801621344;-31.0374084890469,-52.4838438162102;28.0944435741547,-58.2750768122666;18.3970126861515,48.5920982378818;-19.5215744832488,-7.33393904114735;-60.8414503911938,15.5539388242290;46.9498522520006,17.7060090883978;-29.3266608393915,18.8874142536604;-28.1094847324876,52.9484769826052;-66.0324580639326,-28.5874993565832;16.7055570530325,-46.8346897562081;-14.9189668598841,-74.2417493795593;42.8552652265820,-39.5859154983543;26.3507145241858,49.4515323957145;-12.0309445868776,-4.52494211064166;-62.0630131973175,7.64775236610042;-36.6318973496899,-24.0061002516836;-4.53735286063732,-34.9048165820909;67.2981552125186,21.5791430963785;-23.0565300842135,-53.0366366230219;-23.8025442720866,24.6739621149354;-59.1603900887728,41.5513117388991;62.5902294640123,-10.1582085645302;16.0033601290985,5.01417710341792;-67.6449233765971,-20.7516869790786;7.83895136177013,72.9935542653020;-0.826382771885093,-27.8175956175389;-21.8305993620089,47.9910987940474;-17.6628453058907,19.5452101610490;70.8169082870940,14.3945500807536;52.7102064758139,41.0737479150583;21.1007545627210,-40.1502123008772;30.5983588896414,56.2307297925022;-68.3295011127666,2.67471212842579];
VNAx = VNA(:,1);
VNAy = VNA(:,2);
SizeXY = size(VNA,1);
plot(VNAx, VNAy, '*', 'color', 'k', 'markersize', 12);
for i = 1:SizeXY
centers_nodeXY = [VNAx(i), VNAy(i)];
end
%%
% Find Branch Roots on the Edge of the Boundary
ini_branch_dist = nan(numel(VNAx));
CenterEdge_threshold = 80;
for i = 1:SizeXY
for j = 1:SizeCenter_Plane
ini_branch_dist(i,j) = sqrt((CenterPlane_x(j)-VNAx(i)).^2 + (CenterPlane_y(j)-VNAy(i)).^2);
end
end
% find the points that are equal to CenterEdge_threshold:
i2keepEdge_Center = find(ini_branch_dist == CenterEdge_threshold);
% put those into one pair of arrays
keep_x = VNAx(i2keepEdge_Center);
keep_y = VNAy(i2keepEdge_Center);
% #AD: show which markers are "keep"
plot(keep_x, keep_y, 'mo', 'MarkerfaceColor', 'm')
% and the others into another pair of arrays
x_remaining = VNAx;
y_remaining = VNAy;
x_remaining(i2keepEdge_Center) = [];
y_remaining(i2keepEdge_Center) = [];
BranchRoots = [keep_x, keep_y];
RemainingVNs = [x_remaining, y_remaining];
%%
% Computes distances between BranchRoots along circle edge and all remaining vein nodes
distRoot_Branch = pdist2(BranchRoots, RemainingVNs);
% #AD: What do the distances look like?
figure()
imagesc(distRoot_Branch)
% axis equal
cb = colorbar();
ylabel(cb, 'Distance')
axes(ax) % make the original ax current again
%%
% Split RemainingVNs into groups that are nearest to Starting Root VNs in BranchRoots
[~, minRowIdx] = min(distRoot_Branch,[],1);
[GroupID, GroupList] = findgroups(minRowIdx);
Root_BranchNeighborGroups = splitapply(@(x){x},RemainingVNs,GroupID(:));
SizeRoot_BranchNG = size(Root_BranchNeighborGroups,1);
% #AD: Color each group for confirmation
cellfun(@(xy,colr) plot(ax, xy(:,1),xy(:,2), 'o', 'Color', colr, 'MarkerFaceColor', colr), ...
Root_BranchNeighborGroups, ...
mat2cell(jet(numel(Root_BranchNeighborGroups)),ones(numel(Root_BranchNeighborGroups),1), 3));
%%
%Create Branch Matrices
Branch1 = [BranchRoots(1,:);Root_BranchNeighborGroups{1}];
Branch1x = Branch1(:,1);
Branch1y = Branch1(:,2);
SizeofBranch1 = size(Branch1,1);
%%
% Find Branch 1 Terminating VNs or Branch Ends
% #AD in the line below, isn't numel(Branch1x) the same as SizeOfBranch1? Why not just use the variable?
Branch1_dist = nan(numel(Branch1x)); % places nans on the diagonal after distance calculation
VN_VN_proximity = 8;
for i = 1:SizeofBranch1 % looping from largest index to avoid calculating the size of the elim_dist2 matrix without pying the price of dynamic growing of a matrix
for j = 1:(i-1) % distance symmetric (l(i1,i2) == l(i2,i1)) so calculate it once
Branch1_dist(i,j) = sqrt((Branch1x(i)-Branch1x(j)).^2 + (Branch1y(i)-Branch1y(j)).^2);
Branch1_dist(j,i) = Branch1_dist(i,j);
end
end
% #AD See which coordinate is which (label with numbers, requires file from file exchange)
% This is just to match the columns/rows of Branch1_dist onto the points, not really important.
% labelpoints(Branch1x,Branch1y,1:numel(Branch1x), 'NE', .1, 'FontSize', 14) % see file exchange
% What does this distance matrix look like?
% NOTE: This is the upper-most group of points (dark blue)
figure()
imagesc(Branch1_dist)
% axis equal
cb = colorbar();
ylabel(cb, 'Distance')
axes(ax) % make the original ax current again
% find the points that have exactly 1 nearest neighbor
% #AD replacing this line with the ones below (you can rename variables, of course)
% i2keepB1_B1 = find(min(Branch1_dist) >= VN_VN_proximity); % This is the line of code I am trying to figure out.
% #AD if you want to find the nodes that have 1 and only 1 neighbor at the exact distance of 8 units,
% Important: due to floating point roundoff error, the value Branch1_dist(2,1) isn't actually 8 as it
% appears in the command window! For example, Branch1_dist(2,1)-8 = 1.4211e-14.
% To get around this, we'll accept any distances that are very close to 8 with less than 0.000001 error.
isCloseTo_VN_VN_proximity = (abs(Branch1_dist - VN_VN_proximity)) <= 0.000001;
% Find nodes that only have 1 neighbor at ~8 units distance
hasOnlyOneNN = sum(isCloseTo_VN_VN_proximity,1)==1;
% Draw red boxes around the selected nodes
plot(ax, Branch1x(hasOnlyOneNN), Branch1y(hasOnlyOneNN), 'rs','markersize',18,'linewidth',2)
% #AD I didn't look any further than this.
Vance Blake on 19 Jun 2020
Hey Adam, thank you for getting back to me and thanks so much for this solution! Also thanks for the step by step enhancements and explanations.
Adam Danz on 19 Jun 2020

### Categories

Find more on Annotations in Help Center and File Exchange

R2019a

### Community Treasure Hunt

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

Start Hunting!