Subtracting XYZ coordinates defined by matrices

I have two matrices, A (96100x3) and B(16416x3). I need to subtract the XYZ coordinates defined by the columns of B from those of A. But as they are of different resolutions, the difference or setxor is still giving sort of union of these two matrices.

Is there any way I can subtract the region bounded by B from A?

TIA

[Image Analyst]

Please give an example with A having 5 rows and B having 3 rows and give the expected/desired output.

[Isaac John]

I am sorry, but I didnt understand your question. My problem is to subtract the XYZ point cloud described by the B matrix from the XYZ point cloud described by matrix A. Now, the coordinates in B closely matches with some of those in A, but is not exact, and hence setxor(A,B) is not executable. Is there a way I can go ahead with the problem?

[Steven Lord]

What exactly would you expect the result to be if I gave you the following points? How did you compute those expected results? Be specific in describing your process.

If these are not valid values for A and B, explain why they are invalid.

A = [6 1 2; 7 2 7; 1 4 7; 7 7 4; 5 7 6]
B = [2 8 1; 4 9 8; 9 6 9]

[Isaac John]

A matrix would be as it is, in this case, A = [6 1 2; 7 2 7; 1 4 7; 7 7 4; 5 7 6];

B matrix: [6.1 1.2 2.3; 1.2 4.1 7.4];

These are how I expect the matrices to be. I intend to subtract B from A. Notice how the first and second row of B almost matches with first and third rows of A. So, after the subtraction, I am expecting the resultant to be C = [7 2 7; 7 7 4; 5 7 6]. But due to the fact that the rows in B are exactly not matching with A due to the decimal places, I am unable to proceed any further.

I beleive this particular example can be solved with round/floor/ciel commands, but I dont think it is as straight forward for the A and B matrices that I have posted in the OP.

[Image Analyst]

Isaac, so you have this:

A = [6 1 2; 7 2 7; 1 4 7; 7 7 4; 5 7 6]
B = [6.1 1.2 2.3; 1.2 4.1 7.4]
C = [7 2 7; 7 7 4; 5 7 6]

A =

6 1 2

7 2 7

1 4 7

7 7 4

5 7 6

B =

6.1 1.2 2.3

1.2 4.1 7.4

C =

7 2 7

7 7 4

5 7 6

Can you explain how you get each of the 3 rows of C? Why would you want 3 rows from C when A had 5 rows and B had 2 rows? C does not even have any fractional parts - it's all integers even though B is not integers. Explain your thought process.

[Isaac John]

Thank you for your patience.

The idea is two subtract the point cloud defined by B from A. Normally, I would have gone for C = setxor(A,B,'rows'), if the rows in B are exactly matching A. Had B been [6 1 2; 1 4 7] instead of the being [6.1 1.2 2.3; 1.2 4.1 7.4], C = C = setxor(A,B,'rows') would have yielded C = [7 2 7; 7 7 4; 5 7 6], which is basically the point cloud B subtracted from A.

But the case that I have is, the given rows of B are not exactly matching with A and hence, setxor(A,B,'rows') does not yield the desired answer.

For instance let A correspond to the point cloud corresponding to a sphere and B correspond to a smaller, concentric sphere. The result (equivalent to setxor) that I am seeking is a hollow sphere, that is a resultant of subracting the smaller sphere from the larger sphere.

Answer 1

[Adam Danz]

I think what you're asking is how to remove rows of A that are nearly equal to rows of B. Removal is different from subtraction and setxor does the former, not the latter.

The question is, what do you define as "near"?

I used pdist2() to compute the pairwise distances between all points in A and B and there is no obvious cluster of distances around 0 (note, pdist2 will be slow due to large variable sizes!).

% Read-in the data
A = readmatrix('A.xlsx');
B = readmatrix('B.xlsx'); 

% Compute pair-wise distance.  
% This will probably take several seconds or more! size: 96100 x 16416
pd = pdist2(A,B);  

% Plot distribution of pair-wise distances
figure()
histogram(pd(:), 20)
xlabel('paired distance')
ylabel('count')
title('All points')

% Show +/- 1 std
mu = mean(pd(:)); 
sd = std(pd(:)); 
xline(mu-sd,'k-')
xline(mu+sd,'k-')

When we look at distances around 0:.2, there is still no clear cluster that would support an objective definition of "near".

figure()
histogram(pd(pd<.2))
xlabel('paired distance')
ylabel('count')
title('Points at distances less than 0.2')

The minimum distance is

min(pd(:))
% ans =
%       0.000316512059369112

If there were a clear definition of "near" you could use that as a threshold to detect and remove "near" coordinates from the pd pair-wise distance matrix. But that isn't clear from the data.

If, on the other hand, you're looking for coordinates that are nearly 1 unit apart in distanced, then the data above do suggest a clear, more-objective way to detect and eliminate paired coordinates at nearly 1-unit distance, though you'll still need to subjectively select a threshold distance from 1.

For example, any pairwise distance within 0.05 of 1 identifies a row of A to be removed.

rmIdx = abs(pd-1) < .05; 
A(any(rmIdx,2)) = []; 

[Isaac John]

Thank you for your valuable help. Its is reall shame that there is not a well defined threshold to define what is near or not. But I think the data in the above matrices are a constant 1 unit apart. Like as you suggested, I still might be able to come up with a defenition of threshold, which I can try.

One more thing, it appears that the pdist2() is not usable and I am getting the "Requested 96100x16416 (11.8GB) array exceeds maximum array size preference..." error. Is there a way over this?

[Image Analyst]

There can be no well defined threshold because it's a judgment call.

Anyway, pdist2() will not tell you if the points of A are inside the outer shell of points B - it only tells you if they're close, which is not the same thing. For example if the points of A were a shell of radius 20 plus one point at the very center, and the points of B were a shell of radius 19, pdist2() will not tell you that the point of A that's at the center is inside point set B.

You'll have to use

DT = delaunayTriangulation(P)

and then use the isInterior method to determine if the points in A are inside the Delaunay triangulation of set B, and if they are, remove them.

[Adam Danz]

The idea is to minimized any subjective decisions by using the data to set thresholds. This requires exploration of the data.

Here's another visualization of your data zoomed into the section of overlap between A (blue o) and B (red x). This view is from the "top" where each marker is actually a column of markers in 3D.

A = readmatrix('A.xlsx');
B = readmatrix('B.xlsx'); 
plot3(A(:,1), A(:,2), A(:,3), 'o')
hold on
plot3(B(:,1), B(:,2),B(:,3), 'rx')
view(2)
axis equal
xlim([-.05 .25])
ylim([.33 .88])

This is only a small portion of your entire dataset.

And with this view we also see that the mean distance of 1 between points in A and B is irrelevant because you're only concerned with the portion of overlap and those distances a lot closer than 1.

What you need to do is repeat the pdist2 calculations but only for the section of overlap. That will greatly reduce the number of data points and will prevent the maximum array size error you got.

1. Get the min & max values for X and Y coordinates in B
2. Isolate the section of A based on those ranges (add a little space)
3. Repeat the pdist2 calculations I showed you and plot the distribution of distances. The distribution should be much closer to 0.

[Isaac John]

Great job Adam, much appreciated.

Answer 2

[Image Analyst]

"For instance let A correspond to the point cloud corresponding to a sphere and B correspond to a smaller, concentric sphere. The result (equivalent to setxor) that I am seeking is a hollow sphere, that is a resultant of subracting the smaller sphere from the larger sphere."

Sounds like you want convhulln(). Use it on one cloud and then test all the points of one to see if the point is inside the convex hull of the other.