Select only those nodes of an rx3 matrix placed at a distance 'd' from a node P (1x3)

Hi! Is there any way to extract the coordinates of the nodes of the similar circle on the left?

I had thought of something related to the distance 'd' between one node and another of the similar circle but there are nodes that may be at a greater distance from the set 'd' (see nodes A and B). So, using this procedure, it would be necessary first to add nodes between A and B for example. Any methods that could be used?

load M
P = [25.9349 -15.0445 77.3427];

figure
plot3(P(:,1),P(:,2),P(:,3),'k.','Markersize',20);
hold on
plot3(M(:,1),M(:,2),M(:,3),'r.','Markersize',10);
hold off
axis equal

Answer 1

[Voss]

load M
figure
plot3(M(:,1),M(:,2),M(:,3),'r.')
view(2)

% keep only points whose x-coordinate is less than 28:
idx = M(:,1) < 28;
M = M(idx,:);

figure
plot3(M(:,1),M(:,2),M(:,3),'r.')

[Alberto Acri]

It is a good solution, however, I would like to select nodes taking into account P and the fact that the nodes to be selected are close to each other and form a geometry similar to a circle/ellipse.

[Voss]

load M
P = [25.9349 -15.0445 77.3427];

d = sum((M-P).^2,2); % (squared) distance to point P

hist(d) % histogram to help pick distance threshold

d_thresh = 6; % distance threshold

idx = d < d_thresh;
M = M(idx,:);

plot3(M(:,1),M(:,2),M(:,3),'r.')
hold on
plot3(P(1),P(2),P(3),'ko')

[Alberto Acri]

Is this ok as a solution though if I have similar data like this?

load M1
P = [28.73645 -15.2786 69.24005];

figure
plot3(P(:,1),P(:,2),P(:,3),'k.','Markersize',20);
hold on
plot3(M1(:,1),M1(:,2),M1(:,3),'r.','Markersize',10);
hold off
axis equal

[Voss]

Use a different distance threshold for this data set.

load M1
P = [28.73645 -15.2786 69.24005];

figure
plot3(P(:,1),P(:,2),P(:,3),'k.','Markersize',20);
hold on
plot3(M1(:,1),M1(:,2),M1(:,3),'r.','Markersize',10);
hold off
axis equal
view(2)

d = sum((M1-P).^2,2); % (squared) distance to point P

hist(d,100) % histogram (w/ 100 bins) to help pick distance threshold

d_thresh = 5; % distance threshold

idx = d < d_thresh;
M1 = M1(idx,:);

plot3(M1(:,1),M1(:,2),M1(:,3),'r.')
hold on
plot3(P(1),P(2),P(3),'ko')
view(2)

You could try to write a general algorithm - maybe find the lowest-valued histogram bin whose count is zero and use that as the (squared) distance threshold - but that would depend on the width of the bins and I don't think it would be a good enough approach to work in every possible case. And, obviously if your data represents two rings/circles that overlap, then the distance-thresholding approach won't work at all in that case. Maybe someone else has a better idea.

[Alberto Acri]

Ok thanks for the solutions. How can I find the empty bin?

I don't have overlapping circles so I could go with the idea you provided.

[Voss]

Actually, it's not only overlapping rings that would be a problem for the method I've shown above. Any time you have a situation where some points of the ring you want to keep are further from P than some points you don't want to keep, it's not going to work. For instance, if the desired ring is elongated away from the undesired ring and those points far out from P are farther than the closest points in the undesired ring are from P, that won't work.

[Voss]

To illustrate the problem with a simple example. Let's say the P and M are these:

P = [28.73645 -15.2786 69];

M = [ ...
    29.5  -15    69; ...
    29.3  -14.8  69; ...
    29    -14.6  69; ...
    28.75 -14.75 69; ...
    28.5  -15    69; ...
    28.25 -15.2  69; ...
    28    -15.5  69; ...
    28.25 -15.65 69; ...
    28.5  -15.7  69; ...
    28.75 -15.6  69; ...
    29    -15.5  69; ...
    29.25 -15.25 69; ...
    28.25 -14.75 69; ...
    28.3  -14.7 69; ...
    28.2  -14.8 69];

And the last three rows of M are the ones you want to discard.

plot3(M(:,1),M(:,2),M(:,3),'r.')
hold on
plot3(M(end-2:end,1),M(end-2:end,2),M(end-2:end,3),'kx')
plot3(P(1),P(2),P(3),'ko')
view(2)

The distance-threshold approach won't work because those three points are closer to P than some points you want to keep

d = sum((M-P).^2,2); % (squared) distance to point P
disp(d)

    0.6606
    0.5466
    0.5300
    0.2796
    0.1335
    0.2428
    0.5914
    0.3746
    0.2335
    0.1035
    0.1185
    0.2646
    0.5161
    0.5253
    0.5168

Notice that the last three elements (the squared distances to the points to discard) are less than the first element (squared distance to a point to keep).

Therefore, there is no distance threshold you can apply that will keep the first point and discard the last three.

You'll potentially have this problem any time the ring is elongated (more like an ellipse with high eccentricity than a circle).

[Alberto Acri]

Okay thanks for the explanation. Actually I should not have very elliptical shapes. If I still wanted to take advantage of your code how could I determine 'd_thresh' automatically without having to check it in hist?

[Voss]

Actually, the rings don't have to be very elliptical for that method not to work; if there is any part of another ring near enough to a 'flat' side of the ring you want, then those points would incorrectly be included as well. (And, even if there is a distance threshold that would work in a given case, I don't know if there's a very reliable way to find it.)

So here's an alternate method that I think will work better - and it's more in line with the approach you initially conceived of. The idea is to find points that are close to each other (not to find points that are close to P, which was the previous approach I suggested).

% (squared) distance threshold used to determine whether two points are
% 'connected', i.e., close enough to each other to be considered part of
% the same ring:
d_thresh = 0.25;

% note that the same d_thresh is used in both cases here, so maybe (?) the
% same value will work in all cases you have

load M1
P = [28.73645 -15.2786 69.24005];
M_original = M1;

% start with the point in M1 closest to P,
[~,min_idx] = min(sum((M1-P).^2,2));

% and find which points in M1 are "connected" to that point:
is_connected = get_connected_points(M1,M1(min_idx,:),d_thresh);

% matrix of connected points:
M_final = M1(is_connected,:);

% plot results:
figure
plot3(P(1),P(2),P(3),'ok')
hold on
plot3(M_original(:,1),M_original(:,2),M_original(:,3),'r.')
plot3(M_final(:,1),M_final(:,2),M_final(:,3),'g.')

load M
P = [25.9349 -15.0445 77.3427];
M_original = M;

% start with the point in M closest to P,
[~,min_idx] = min(sum((M-P).^2,2));

% and find which points in M are "connected" to that point:
is_connected = get_connected_points(M,M(min_idx,:),d_thresh);

% matrix of connected points:
M_final = M(is_connected,:);

% plot results:
figure
plot3(P(1),P(2),P(3),'ok')
hold on
plot3(M_original(:,1),M_original(:,2),M_original(:,3),'r.')
plot3(M_final(:,1),M_final(:,2),M_final(:,3),'g.')

function is_connected = get_connected_points(M,P,d_thresh,is_connected)

% inputs:
%            M: Nx3 matrix of points
%            P: reference point
%       thresh: squared distance threshold
% is_connected: Nx1 logical vector saying whether each point has been found
%               to be connected (i.e., closer than thresh) to some point in M so far;
%               this is not given initially but is used in recursive calls
% output:
% is_connected: (updated) is_connected vector

if nargin < 4
    % if no is_connected vector is given, initialize to all false
    % (no points have been found to be connected yet at the beginning)
    is_connected = false(size(M,1),1);
elseif all(is_connected)
    % if all points have been found to be connected, nothing left to do
    return
end

% d: squared distance between each point in M and reference point P
% d is Inf for points already found to be connected so that they are not "found" again
d = Inf(size(M,1),1);

% for those points still in play (i.e., not found to be connected to any
% other point (yet)), calculate their squared distance from reference point P
d(~is_connected) = sum((M(~is_connected,:)-P).^2,2);

% find the indices of points closer than d_thresh to the reference point P
idx = find(d < d_thresh);

% mark those points as connected
is_connected(idx) = true;

% for each of those newly-found connected points,
for ii = 1:numel(idx)
    % find the points connected to them, i.e., run this function again,
    % with each newly-found point as the reference point P, updating and
    % returning the is_connected vector each time
    is_connected = get_connected_points(M,M(idx(ii),:),d_thresh,is_connected);
end
end

[Alberto Acri]

Ok thanks for the help @Voss! This solution might be the most useful.

[Voss]

You're welcome!