Fast method to remove redundant points in Matlab

I merge two pointCloud objects in Matlab, let's say pc1 and pc2. pc1 is the reference cloud, i.e. all points from pc2 which are equal or very close to points from pc1 need to be removed before the clouds are combined.

Clarifications:

I know the function pcmerge, which is almost doing what I want - but I definitely need to remove redundant points and averaging those points is not an option
The point cloud size is about 500,000 each and I have to compare many (100) of them. That's why the speed is important.
I would prefer to be able to define a radius around each point of pc1 to give a criterion for "being redundant". But in favor of the speed, some simplifications are okay (cf. my second solution approach).

Solution approaches:

A working but really slow solution is to look for each point in pc2 for its nearest neighbor:

function [ pc ] = pcaddcloud( pc1, pc2, res )

limits = overlapRange(pc2, pc1);
pc1idx = findPointsInROI(pc2, limits);
pc2Overlap = select(pc2, pc1idx);
idx = findPointsInROI(pc1, limits);
pc1Overlap = select(pc1, idx);
endi = pc2Overlap.Count;
pc2Overlap = pc2Overlap.Location;
for i=1:endi
    [idx, ~] = findNeighborsInRadius(pc1Overlap, pc2Overlap(i,:), res);
    % keep only indices of redundant points to delete them later
    if isempty(idx)
        pc1idx(i) = 0;
    end
end
pc1idx(pc1idx==0) = [];
pc2 = pc2.Location;
pc2(pc1idx,:) = [];
pc = pointCloud([pc1.Location; pc2]);
end

% Compute the bounding box of overlapped region (from pcmerge)
function rangeLimits = overlapRange(pcA, pcB)
xlimA = pcA.XLimits;
ylimA = pcA.YLimits;
zlimA = pcA.ZLimits;

xlimB = pcB.XLimits;
ylimB = pcB.YLimits;
zlimB = pcB.ZLimits;

if (xlimA(1) > xlimB(2) || xlimA(2) < xlimB(1) || ...
        ylimA(1) > ylimB(2) || ylimA(2) < ylimB(1) || ...
        zlimA(1) > zlimB(2) || zlimA(2) < zlimB(1))
    % No overlap
    rangeLimits = [];
else
    rangeLimits = [ min(xlimA(1),xlimB(1)), max(xlimA(2),xlimB(2)); ...
        min(ylimA(1),ylimB(1)), max(ylimA(2),ylimB(2)); ...
        min(zlimA(1),zlimB(1)), max(zlimA(2),zlimB(2))];
end
end

I have a faster solution (still slow, but faster than solution 1) working with alpha shapes: I define a hull around pc1 and decide if the points of pc2 are inside or not. Disadvantage: points, that are only "slightly outside" (i.e. close to points of pc1 but not inside the alpha shape) are not detected as redundant.

function [ pc ] = pcaddcloud( pc1, pc2 )

limits = overlapRange(pc2, pc1);
pc2 = pc2.Location;
pc1 = pc1.Location;
%seems to be faster than findPointsInROI:
pc2Overlap = pc2(pc2(:,1)>=limits(1,1)&pc2(:,1)<=limits(1,2) ...
    &pc2(:,2)>=limits(2,1)&pc2(:,2)<=limits(2,2)...
    &pc2(:,3)>=limits(3,1)&pc2(:,3)<=limits(3,2),:);
pc2idx = find(pc2(:,1)>=limits(1,1)&pc2(:,1)<=limits(1,2) ...
    &pc2(:,2)>=limits(2,1)&pc2(:,2)<=limits(2,2)...
    &pc2(:,3)>=limits(3,1)&pc2(:,3)<=limits(3,2));
pc1Overlap = pc1(pc1(:,1)>=limits(1,1)&pc1(:,1)<=limits(1,2) ...
    &pc1(:,2)>=limits(2,1)&pc1(:,2)<=limits(2,2)...
    &pc1(:,3)>=limits(3,1)&pc1(:,3)<=limits(3,2),:);

shape = alphaShape(double(pc1Overlap));
in = inShape(shape, double(pc2Overlap));
pc2idx(~in) = [];
pc2(pc2idx,:) = [];
pc = pointCloud([pc1; pc2]);

end

% Compute the bounding box of overlapped region (from pcmerge)
function rangeLimits = overlapRange(pcA, pcB)
xlimA = pcA.XLimits;
ylimA = pcA.YLimits;
zlimA = pcA.ZLimits;

xlimB = pcB.XLimits;
ylimB = pcB.YLimits;
zlimB = pcB.ZLimits;

if (xlimA(1) > xlimB(2) || xlimA(2) < xlimB(1) || ...
        ylimA(1) > ylimB(2) || ylimA(2) < ylimB(1) || ...
        zlimA(1) > zlimB(2) || zlimA(2) < zlimB(1))
    % No overlap
    rangeLimits = [];
else
    rangeLimits = [ min(xlimA(1),xlimB(1)), max(xlimA(2),xlimB(2)); ...
        min(ylimA(1),ylimB(1)), max(ylimA(2),ylimB(2)); ...
        min(zlimA(1),zlimB(1)), max(zlimA(2),zlimB(2))];
end
end

I'm looking forward to your ideas! Please do not hesitate to ask for more information, if needed - I am new to this platform. Thank you!

Solution

You can use ismembertol with ByRows option to detect redundant points. But consider that instead of spherical neighborhood it uses cubic neighborhood. Assuming that you have two matrices pc1 ,pc2 each one has 3 columns and a tolerance tol:

idx = ismembertol(pc2, pc1, tol,'ByRows', true, 'DataScale' , 1);
result = [pc1; pc2(~idx,:)];