planePosition

PLANEPOSITION Compute position of a point on a plane.

   COORDS2D = planePosition(POINT, PLANE)
   POINT has format [X Y Z], and plane has format
   [X0 Y0 Z0  DX1 DY1 DZ1  DX2 DY2 DZ2], where :
   - (X0, Y0, Z0) is a point belonging to the plane
   - (DX1, DY1, DZ1) is a first direction vector
   - (DX2, DY2, DZ2) is a second direction vector
   Result COORDS2D has the form [XP YP], with XP and YP the coordinates of
   the point in the coordinate system of the plane.

   [COORDS2D, Z] = planePosition(POINT, PLANE)
   Also returns the coordinates along the normal of the plane. When the
   point is within the plane, this coordinate should be zero.

   Example
     plane = [10 20 30  1 0 0  0 1 0];
     point = [13 24 35];
     pos = planePosition(point, plane)
     pos = 
         3   4

   Example
     % plane with non unit direction vectors
     p0 = [30 20 10]; v1 = [2 1 0]; v2 = [-2 4 0];
     plane = [p0 v1 v2];
     pts = [p0 ; p0 + v1 ; p0 + v2 ; p0 + 3 * v1 + 2 * v2];
     pos = planePosition(pts, plane)
     pos = 
          0     0
          1     0
          0     1
          3     2


   See also
     geom3d, planes3d, points3d, planePoint

0001 function [pos, zp] = planePosition(point, plane)
0002 %PLANEPOSITION Compute position of a point on a plane.
0003 %
0004 %   COORDS2D = planePosition(POINT, PLANE)
0005 %   POINT has format [X Y Z], and plane has format
0006 %   [X0 Y0 Z0  DX1 DY1 DZ1  DX2 DY2 DZ2], where :
0007 %   - (X0, Y0, Z0) is a point belonging to the plane
0008 %   - (DX1, DY1, DZ1) is a first direction vector
0009 %   - (DX2, DY2, DZ2) is a second direction vector
0010 %   Result COORDS2D has the form [XP YP], with XP and YP the coordinates of
0011 %   the point in the coordinate system of the plane.
0012 %
0013 %   [COORDS2D, Z] = planePosition(POINT, PLANE)
0014 %   Also returns the coordinates along the normal of the plane. When the
0015 %   point is within the plane, this coordinate should be zero.
0016 %
0017 %   Example
0018 %     plane = [10 20 30  1 0 0  0 1 0];
0019 %     point = [13 24 35];
0020 %     pos = planePosition(point, plane)
0021 %     pos =
0022 %         3   4
0023 %
0024 %   Example
0025 %     % plane with non unit direction vectors
0026 %     p0 = [30 20 10]; v1 = [2 1 0]; v2 = [-2 4 0];
0027 %     plane = [p0 v1 v2];
0028 %     pts = [p0 ; p0 + v1 ; p0 + v2 ; p0 + 3 * v1 + 2 * v2];
0029 %     pos = planePosition(pts, plane)
0030 %     pos =
0031 %          0     0
0032 %          1     0
0033 %          0     1
0034 %          3     2
0035 %
0036 %
0037 %   See also
0038 %     geom3d, planes3d, points3d, planePoint
0039 %
0040 
0041 % ------
0042 % Author: David Legland
0043 % E-mail: david.legland@inrae.fr
0044 % Created: 2005-02-21
0045 % Copyright 2005-2024 INRA - TPV URPOI - BIA IMASTE
0046 
0047 % size of input arguments
0048 npl = size(plane, 1);
0049 npt = size(point, 1);
0050 
0051 % check inputs have compatible sizes
0052 if npl ~= npt && npl > 1 && npt > 1
0053     error('geom3d:planePosition:inputSize', ...
0054         'plane and point should have same size, or one of them must have 1 row');
0055 end
0056 
0057 % origin and direction vectors of the plane(s)
0058 p0 = plane(:, 1:3);
0059 v1 = plane(:, 4:6);
0060 v2 = plane(:, 7:9);
0061 
0062 % Principle
0063 %
0064 % for each (recentered) point, we need to solve the system:
0065 %   s * v1x + t * v2x + u * v3x  = px
0066 %   s * v1y + t * v2y + u * v3y  = py
0067 %   s * v1z + t * v2z + u * v3z  = pz
0068 % Assuming the point belong to the place, the value of u is 0. The last
0069 % equation is kept only for homogeneity.
0070 % We rewrite in matrix form:
0071 %   [ v1x v2x v3x ]   [ s ]   [ xp ]
0072 %   [ v1y v2y v3y ] * [ t ] = [ yp ]
0073 %   [ v1z v2z v3z ]   [ u ]   [ zp ]
0074 % Or:
0075 %   A * X = B
0076 % We need to solve X = inv(A) * B. In practice, we use the b/A notation,
0077 % obtained after using transpositon on the A\b notation.
0078 %   X' = (inv(mat) * bsxfun(@minus, point, p0)')';
0079 %   X' = bsxfun(@minus, point, p0) * inv(mat');
0080 %   X' = bsxfun(@minus, point, p0) / mat';
0081 %   (and we compute the transpose of the basis transform)
0082 
0083 % Compute dot products with direction vectors of the plane
0084 if npl == 1
0085     % we have npl == 1 and npt > 1
0086     % build transpose of matrix that changes plane basis to global basis
0087     mat = [v1 ; v2 ; cross(v1, v2, 2)];
0088     tmp = bsxfun(@minus, point, p0) / mat;
0089     pos = tmp(:, 1:2);
0090     zp = tmp(:,3);
0091 
0092 else
0093     % NPL > 0 -> iterate over planes
0094     % Number of points can be either 1, or the same number as planes
0095     pos = zeros(npl, 2);
0096     zp = zeros(npl, 1);
0097     for ipl = 1:npl
0098         mat = [v1(ipl,:) ; v2(ipl,:) ; cross(v1(ipl,:), v2(ipl,:), 2)];
0099         % choose either point with same index, or the single point
0100         ind = min(ipl, npt);
0101         tmp = bsxfun(@minus, point(ind,:), p0(ipl,:)) / mat;
0102         pos(ipl,:) = tmp(1,1:2);
0103         zp(ipl) = tmp(1,3);
0104     end
0105 end
0106 
0107 % % old version (:
0108 % s = dot(point-p0, d1, 2) ./ vectorNorm3d(d1);
0109 % t = dot(point-p0, d2, 2) ./ vectorNorm3d(d2);