projPointOnPlane

PROJPOINTONPLANE Return the orthogonal projection of a point on a plane.

   PT2 = projPointOnPlane(PT1, PLANE);
   Compute the (orthogonal) projection of point PT1 onto the plane PLANE,
   given as [X0 Y0 Z0  VX1 VY1 VZ1  VX2 VY2 VZ2] (origin point, first
   direction vector, second directionvector).
   
   The function is fully vectorized, in that multiple points may be
   projected onto multiple planes in a single call, returning multiple
   points. With the exception of the second dimension (where
   SIZE(PT1,2)==3, and SIZE(PLANE,2)==9), each dimension of PT1 and PLANE
   must either be equal or one, similar to the requirements of BSXFUN. In
   basic usage, point PT1 is a [N*3] array, and PLANE is a [N*9] array
   (see createPlane for details). Result PT2 is a [N*3] array, containing
   coordinates of orthogonal projections of PT1 onto planes PLANE. In
   vectorised usage, PT1 is an [N*3*M*P...] matrix, and PLANE is an
   [X*9*Y...] matrix, where (N,X), (M,Y), etc, are either equal pairs, or
   one of the two is one.

   See also:
   planes3d, points3d, planePosition, intersectLinePlane

0001 function point = projPointOnPlane(point, plane)
0002 %PROJPOINTONPLANE Return the orthogonal projection of a point on a plane.
0003 %
0004 %   PT2 = projPointOnPlane(PT1, PLANE);
0005 %   Compute the (orthogonal) projection of point PT1 onto the plane PLANE,
0006 %   given as [X0 Y0 Z0  VX1 VY1 VZ1  VX2 VY2 VZ2] (origin point, first
0007 %   direction vector, second directionvector).
0008 %
0009 %   The function is fully vectorized, in that multiple points may be
0010 %   projected onto multiple planes in a single call, returning multiple
0011 %   points. With the exception of the second dimension (where
0012 %   SIZE(PT1,2)==3, and SIZE(PLANE,2)==9), each dimension of PT1 and PLANE
0013 %   must either be equal or one, similar to the requirements of BSXFUN. In
0014 %   basic usage, point PT1 is a [N*3] array, and PLANE is a [N*9] array
0015 %   (see createPlane for details). Result PT2 is a [N*3] array, containing
0016 %   coordinates of orthogonal projections of PT1 onto planes PLANE. In
0017 %   vectorised usage, PT1 is an [N*3*M*P...] matrix, and PLANE is an
0018 %   [X*9*Y...] matrix, where (N,X), (M,Y), etc, are either equal pairs, or
0019 %   one of the two is one.
0020 %
0021 %   See also:
0022 %   planes3d, points3d, planePosition, intersectLinePlane
0023 
0024 %   ---------
0025 %   author : David Legland
0026 %   INRA - TPV URPOI - BIA IMASTE
0027 %   created the 18/02/2005.
0028 %
0029 
0030 %   HISTORY
0031 %   21/08/2006: debug support for multiple points or planes
0032 %   22/04/2013: uses bsxfun for mult. pts/planes in all dimensions (Sven H)
0033 
0034 % Unpack the planes into origins and normals, keeping original shape
0035 plSize = size(plane);
0036 plSize(2) = 3;
0037 [origins, normals] = deal(zeros(plSize));
0038 origins(:) = plane(:,1:3,:);
0039 normals(:) = crossProduct3d(plane(:,4:6,:), plane(:, 7:9,:));
0040 
0041 % difference between origins of plane and point
0042 dp = bsxfun(@minus, origins, point);
0043 
0044 % relative position of point on normal's line
0045 t = bsxfun(@rdivide, sum(bsxfun(@times,normals,dp),2), sum(normals.^2,2));
0046 
0047 % add relative difference to project point back to plane
0048 point = bsxfun(@plus, point, bsxfun(@times, t, normals));