Selecting points which belong to parallel planes

Hi and thanks for looking

I have a cloud of points, which very roughly form an irregular cylinder. I need to draw a number of parallel planes through it and calculate which points belong to which plane. Then, I need to calculate the squared sum of distances from those points to the principal axis. So, for each plane there will be 1 value (sq sum of distances). I have calculated its principal axis and am not sure how to proceed next. I know that the principal axis is the normal to the planes

Thank you for your help

Answer 1

[Roger Stafford]

It's a matter of simple vector analysis. Let P0 = [x0,y0,z0] be some base point on your principal axis and let n = [nx,ny,nz] be a vector pointing along this axis from P0. Now let P1 = [x1,y1,z1] be some point in your cloud. You wish to find the point P2 = [x2,y2,z2] at which the principal axis intersects the plane which is orthogonal to the axis and contains P1, and you want to find the squared distance between P1 and P2. Then we know that

1) P2-P1 is orthogonal to n

and

2) P2-P0 is equal to some scalar multiple of n

Combining these two facts we can arrive at a solution for P2:

P2 = P0 + dot(P1-P0,n)/dot(n,n)*n;

Then the squared distance between P1 and P2 would be dot(P1-P2,P1-P2).

I should add that it is easy to vectorize the expressions for all the P2 points and the squared distances, making use of the bsxfun function.

[Roger Stafford]

I decided to show you how your problem could be vectorized. Let P0 be a 3 x 1 column vector of the three coordinates of some point on your axis. Let n be a 3 x 1 column vector of a vector pointing from P0 along the axis. Let P1 be a 3 x m array in which each column is one of m points in the cloud. Then do this:

   % Set up the cloud and axis:
   P0 = randn(3,1); % Some point on the axis
   n = randn(3,1); % Vector pointing along the axis from P0
   m = 1000; % The number of points in the cloud
   P1 = randn(3,m); % Each column is a point in the cloud

   % Computation begins:
   P10 = bsxfun(@minus,P1,P0);
   ne  = repmat(n/norm(n),1,m); % Normalized vector, repeated in m columns
   P20 = bsxfun(@times,sum(P10.*ne,1),ne); % The P2-P0 vectors
   P2  = bsxfun(@plus,P0,P20); % The projections of cloud points onto axis
   D2  = sum((P10-P20).^2,1);  % The squared distances of cloud points from projections

The P2 columns will be the projections of your cloud points onto the axis and D2 will be the squares of the distances between the cloud points and their projections.

Answer 2

[Image Analyst]

First of all, rotate the array so that your principal axis is along the Z direction. Then get the center of the rotated cloud: xCenter, yCenter.

Then, assuming you have 3 vectors for x, y, and z, then you can find all points in a slice within "tolerance" of some z value, zSlice, like this (untested):

tolerance = 0.4; % Whatever thickness you want.
indexesInSlice = abs(z - zSlice) < tolerance;
% Then extract x and y that "belong" to that slice
inSliceX = x(indexesInSlice);
inSliceY = y(indexesInSlice);

Now compute "the squared sum of distances from those points to the principal axis."

sumOfDistances = sum(sqrt((inSliceX - xCenter).^2 + (inSliceY - yCenter).^2));  
squaredSumOfDistances = sumOfDistances .^ 2;

Or if you really meant "the sum of squared distances from those points to the principal axis."

sumOfSquaredDistances = sum((inSliceX - xCenter).^2 + (inSliceY - yCenter).^2);

[Tim Navr]

Rotating it to align the principal axis with anything else makes further calculations wrong though. Unfortunately, this solution is no good. Thank you for your time though

[Image Analyst]

I don't see why you say that, but whatever....good luck.

[Tim Navr]

This is only a part of the code and rotation affects the later parts. It worked well though. Thank you