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

3 views (last 30 days)
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;
radii_plane = 80;
radii_vein = 4;
radii_inf_range = 12;
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);
viscircles(center_plane, radii_plane, 'color', 'b');
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)];
viscircles(centers_nodeXY, radii_vein, 'color', 'k');
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];
%%
  1 Comment
Adam Danz
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);

Sign in to comment.

Accepted Answer

Adam Danz
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;
radii_plane = 80;
radii_vein = 4;
radii_inf_range = 12;
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);
viscircles(center_plane, radii_plane, 'color', 'b');
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)];
viscircles(centers_nodeXY, radii_vein, 'color', 'k');
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
% #AD
% 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.

More Answers (0)

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!