projPointOnCylinder

PROJPOINTONCYLINDER Project a 3D point onto a cylinder.

   PTPROJ = projPointOnCircle3d(PT, CYLINDER).
   Computes the projection of 3D point PT onto the CYLINDER. 
   
   Point PT is a 1-by-3 array, and CYLINDER is a 1-by-7 array.
   Result PTPROJ is a 1-by-3 array, containing the coordinates of the
   projection of PT onto the CYLINDER.

   PTPROJ = projPointOnCircle3d(..., OPT)
   with OPT = 'open' (0) (default) or 'closed' (1), specify if the bases 
   of the cylinder should be included.

   Example
       demoProjPointOnCylinder

   See also
       projPointOnLine3d, projPointOnPlane, projPointOnCircle3d

0001 function ptProj = projPointOnCylinder(pt, cyl, varargin)
0002 %PROJPOINTONCYLINDER Project a 3D point onto a cylinder.
0003 %
0004 %   PTPROJ = projPointOnCircle3d(PT, CYLINDER).
0005 %   Computes the projection of 3D point PT onto the CYLINDER.
0006 %
0007 %   Point PT is a 1-by-3 array, and CYLINDER is a 1-by-7 array.
0008 %   Result PTPROJ is a 1-by-3 array, containing the coordinates of the
0009 %   projection of PT onto the CYLINDER.
0010 %
0011 %   PTPROJ = projPointOnCircle3d(..., OPT)
0012 %   with OPT = 'open' (0) (default) or 'closed' (1), specify if the bases
0013 %   of the cylinder should be included.
0014 %
0015 %   Example
0016 %       demoProjPointOnCylinder
0017 %
0018 %   See also
0019 %       projPointOnLine3d, projPointOnPlane, projPointOnCircle3d
0020 %
0021 
0022 % ---------
0023 % Author: oqilipo
0024 % Created: 2021-04-17, using R2020b
0025 % Copyright 2021
0026 
0027 parser = inputParser;
0028 addRequired(parser, 'pt', @(x) validateattributes(x, {'numeric'},...
0029     {'size',[1 3],'real','finite','nonnan'}));
0030 addRequired(parser, 'cyl', @(x) validateattributes(x, {'numeric'},...
0031     {'size',[1 7],'real','finite','nonnan'}));
0032 capParValidFunc = @(x) (islogical(x) ...
0033     || isequal(x,1) || isequal(x,0) || any(validatestring(x, {'open','closed'})));
0034 addOptional(parser,'cap','open', capParValidFunc);
0035 parse(parser,pt,cyl,varargin{:});
0036 pt = parser.Results.pt;
0037 cyl = parser.Results.cyl;
0038 cap = lower(parser.Results.cap(1));
0039 
0040 % Radius of the cylinder
0041 cylRadius = cyl(7);
0042 % Height of the cylinder
0043 cylBottom = -Inf;
0044 cylHeight = Inf;
0045 if cap == 'c' || cap == 1
0046     cylBottom = 0;
0047     cylHeight = distancePoints3d(cyl(1:3),cyl(4:6));
0048 end
0049 % Create a transformation for the point into a local cylinder coordinate
0050 % system. Align the cylinder axis with the z axis and translate the
0051 % starting point of the cylinder to the origin.
0052 TFM = createRotationVector3d(cyl(4:6)-cyl(1:3), [0 0 1])*createTranslation3d(-cyl(1:3));
0053 % cylTfm = [transformPoint3d(cyl(1:3), TFM) transformPoint3d(cyl(4:6), TFM) cylRadius];
0054 % cylTfm2 = [0 0 0 0 0 cylHeight, cylRadius];
0055 % assert(ismembertol(cylTfm,cylTfm2,'byRows',1,'DataScale',1e1))
0056 
0057 % Transform the point.
0058 ptTfm = transformPoint3d(pt,TFM);
0059 % Convert the transformed point to cylindrical coordinates.
0060 [ptTheta, ptRadius, ptHeight] = cart2cyl(ptTfm);
0061 
0062 if ptRadius <= cylRadius && (ptHeight <= cylBottom || ptHeight >= cylHeight)
0063     % If point is inside the radius of the cylinder but outside its height
0064     if ptHeight <= cylBottom
0065         ptProj_cyl = [ptTheta, ptRadius, 0];
0066     else
0067         ptProj_cyl = [ptTheta, ptRadius, cylHeight];
0068     end
0069 elseif ptRadius > cylRadius && (ptHeight <= cylBottom || ptHeight >= cylHeight)
0070     % If point is outside the cylinder's radius and height
0071     if ptHeight <= cylBottom
0072         ptProj_cyl = [ptTheta, cylRadius, 0];
0073     else
0074         ptProj_cyl = [ptTheta, cylRadius, cylHeight];
0075     end
0076 elseif ptRadius < cylRadius && (ptHeight > cylBottom && ptHeight < cylHeight)
0077     % If point is inside the cylinder's radius and height
0078     deltaRadius = cylRadius - ptRadius;
0079     deltaHeight = cylHeight - ptHeight;
0080     if (deltaRadius < ptHeight && deltaRadius < deltaHeight) || isinf(cylBottom)
0081         % If the distance to the cylinder's surface is smaller than the
0082         % distance to the top and bottom surfaces.
0083         ptProj_cyl = [ptTheta, cylRadius, ptHeight];
0084     else
0085         if ptHeight < deltaHeight
0086             ptProj_cyl = [ptTheta, ptRadius, 0];
0087         else
0088             ptProj_cyl = [ptTheta, ptRadius, cylHeight];
0089         end
0090     end
0091 elseif ptRadius >= cylRadius && (ptHeight > cylBottom && ptHeight < cylHeight)
0092     % If point is outside the radius of the cylinder and inside its height
0093     ptProj_cyl = [ptTheta, cylRadius, ptHeight];
0094 end
0095 
0096 % Convert the projected point back to Cartesian coordinates
0097 ptProj_cart = cyl2cart(ptProj_cyl);
0098 % Transform the projected point back to the global coordinate system
0099 ptProj = transformPoint3d(ptProj_cart,inv(TFM));
0100 
0101 end