geodesicCylinder

GEODESICCYLINDER Computes the geodesic between two points on a cylinder.

   [GEO, GEOLENGTH] = geodesicCylinder(PTS, CYL) computes the geodesic 
   between the two points PTS projected onto the infinite cylinder CYL.  
   PTS is a 2-by-3 array, and CYL is a 1-by-7 array. Result is the 
   polyline GEO (500-by-3 array) [500 = default] containing the  
   coordinates of the geodesic between two projected points. GEOLENGTH 
   contains the analytical length of the geodesic.

   [~, ~, CONGEO, CONGEOLENGTH] = geodesicCylinder(PTS, CYL) provides the
   conjugate geodesic and its analytical length.

   ... = geodesicCylinder(PTS, CYL, 'n', N) defines the number of points
   representing the geodesic and conjugate geodesic.

   Example
       demoGeodesicCylinder

   See also
     drawCylinder, projPointOnCylinder

   Source
     Based on the script 'geodesic.m' by Lei Wang
     https://mathworks.com/matlabcentral/fileexchange/6522

0001 function [geo, geoLength, conGeo, conGeoLength] = geodesicCylinder(pts, cyl, varargin)
0002 %GEODESICCYLINDER Computes the geodesic between two points on a cylinder.
0003 %
0004 %   [GEO, GEOLENGTH] = geodesicCylinder(PTS, CYL) computes the geodesic
0005 %   between the two points PTS projected onto the infinite cylinder CYL.
0006 %   PTS is a 2-by-3 array, and CYL is a 1-by-7 array. Result is the
0007 %   polyline GEO (500-by-3 array) [500 = default] containing the
0008 %   coordinates of the geodesic between two projected points. GEOLENGTH
0009 %   contains the analytical length of the geodesic.
0010 %
0011 %   [~, ~, CONGEO, CONGEOLENGTH] = geodesicCylinder(PTS, CYL) provides the
0012 %   conjugate geodesic and its analytical length.
0013 %
0014 %   ... = geodesicCylinder(PTS, CYL, 'n', N) defines the number of points
0015 %   representing the geodesic and conjugate geodesic.
0016 %
0017 %   Example
0018 %       demoGeodesicCylinder
0019 %
0020 %   See also
0021 %     drawCylinder, projPointOnCylinder
0022 %
0023 %   Source
0024 %     Based on the script 'geodesic.m' by Lei Wang
0025 %     https://mathworks.com/matlabcentral/fileexchange/6522
0026 %
0027 
0028 % ------
0029 % Author: oqilipo
0030 % E-mail: N/A
0031 % Created: 2021-04-17, using MATLAB R2022b
0032 % Copyright 2021-2024
0033 
0034 parser = inputParser;
0035 addRequired(parser, 'pts', @(x) validateattributes(x, {'numeric'},...
0036     {'size',[2 3],'real','finite','nonnan'}));
0037 addRequired(parser, 'cyl', @(x) validateattributes(x, {'numeric'},...
0038     {'size',[1 7],'real','finite','nonnan'}));
0039 addParameter(parser,'n',500, @(x) validateattributes(x, {'numeric'},...
0040     {'scalar','>', 2,'<=', 1e5}));
0041 parse(parser,pts,cyl,varargin{:});
0042 pts = parser.Results.pts;
0043 cyl = parser.Results.cyl;
0044 n = parser.Results.n;
0045 
0046 % Radius of the cylinder
0047 cylRadius = cyl(7);
0048 
0049 % Project points onto the open (infinite) cylinder
0050 ptProj(1,:) = projPointOnCylinder(pts(1,:), cyl, 'open');
0051 ptProj(2,:) = projPointOnCylinder(pts(2,:), cyl, 'open');
0052 
0053 % Create a transformation for the points into the local cylinder coordinate
0054 % system. Align the cylinder axis with the z axis and translate the
0055 % starting point of the cylinder to the origin.
0056 TFM = createRotationVector3d(cyl(4:6)-cyl(1:3), [0 0 1])*createTranslation3d(-cyl(1:3));
0057 % Transform the points.
0058 ptTfm = transformPoint3d(ptProj, TFM);
0059 % Convert the transformed points to cylindrical coordinates.
0060 [ptsTheta, ptsRadius, ptsHeight] = cart2cyl(ptTfm);
0061 assert(ismembertol(ptsRadius(1),ptsRadius(2)))
0062 assert(ismembertol(ptsRadius(1),cylRadius))
0063 
0064 % Copy thetas for the conjugate geodesic
0065 ptsTheta(:,:,2) = ptsTheta;
0066 ptsTheta(1,1,2) = ptsTheta(1,1,2) + 2*pi;
0067 
0068 geoCyl = nan(n,3,size(ptsTheta,3));
0069 arcLength = nan(1,size(ptsTheta,3));
0070 for t = 1:size(ptsTheta,3)
0071     [geoCyl(:,:,t), arcLength(t)] = geoCurve(ptsTheta(:,:,t), cylRadius, ptsHeight, n);
0072 end
0073 
0074 % Select the shortest geodesic
0075 if arcLength(1) <= arcLength(2)
0076     % Transform the geodesics back to the global coordinate system
0077     geo = transformPoint3d(geoCyl(:,:,1), inv(TFM));
0078     conGeo = transformPoint3d(geoCyl(:,:,2), inv(TFM));
0079     geoLength = arcLength(1);
0080     conGeoLength = arcLength(2);
0081 else
0082     % Transform the geodesics back to the global coordinate system
0083     geo = transformPoint3d(geoCyl(:,:,2), inv(TFM));
0084     conGeo = transformPoint3d(geoCyl(:,:,1), inv(TFM));
0085     geoLength = arcLength(2);
0086     conGeoLength = arcLength(1);
0087 end
0088 
0089 end
0090 
0091 function [geo, arcLength] = geoCurve(theta, r, z, n)
0092 % Parametric expression of the geodesic curve
0093 u = linspace(theta(1),theta(2),n)';
0094 geo(:,1) = r*cos(u);
0095 geo(:,2) = r*sin(u);
0096 geo(:,3) = (z(2)-z(1))/(theta(2)-theta(1))*u + (z(1)*theta(2)-z(2)*theta(1))/(theta(2)-theta(1));
0097 if all(isnan(geo(:,3)))
0098     geo(:,3) = linspace(z(1),z(2),n)';
0099 end
0100 arcLength = sqrt(r^2*(theta(2)-theta(1))^2+(z(2)-z(1))^2);
0101 end