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 % E-mail: david.legland@inrae.fr
0027 % Created: 2005-02-18
0028 % Copyright 2005-2024 INRA - TPV URPOI - BIA IMASTE
0029 
0030 % Unpack the planes into origins and normals, keeping original shape
0031 plSize = size(plane);
0032 plSize(2) = 3;
0033 [origins, normals] = deal(zeros(plSize));
0034 origins(:) = plane(:,1:3,:);
0035 normals(:) = crossProduct3d(plane(:,4:6,:), plane(:, 7:9,:));
0036 
0037 % difference between origins of plane and point
0038 dp = bsxfun(@minus, origins, point);
0039 
0040 % relative position of point on normal's line
0041 t = bsxfun(@rdivide, sum(bsxfun(@times,normals,dp),2), sum(normals.^2,2));
0042 
0043 % add relative difference to project point back to plane
0044 point = bsxfun(@plus, point, bsxfun(@times, t, normals));