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 % Created: 2021-04-17, using R2020b
0031 % Copyright 2021
0032 
0033 parser = inputParser;
0034 addRequired(parser, 'pts', @(x) validateattributes(x, {'numeric'},...
0035     {'size',[2 3],'real','finite','nonnan'}));
0036 addRequired(parser, 'cyl', @(x) validateattributes(x, {'numeric'},...
0037     {'size',[1 7],'real','finite','nonnan'}));
0038 addParameter(parser,'n',500, @(x) validateattributes(x, {'numeric'},...
0039     {'scalar','>', 2,'<=', 1e5}));
0040 parse(parser,pts,cyl,varargin{:});
0041 pts = parser.Results.pts;
0042 cyl = parser.Results.cyl;
0043 n = parser.Results.n;
0044 
0045 % Radius of the cylinder
0046 cylRadius = cyl(7);
0047 
0048 % Project points onto the open (infinite) cylinder
0049 ptProj(1,:) = projPointOnCylinder(pts(1,:), cyl, 'open');
0050 ptProj(2,:) = projPointOnCylinder(pts(2,:), cyl, 'open');
0051 
0052 % Create a transformation for the points into the local cylinder coordinate
0053 % system. Align the cylinder axis with the z axis and translate the
0054 % starting point of the cylinder to the origin.
0055 TFM = createRotationVector3d(cyl(4:6)-cyl(1:3), [0 0 1])*createTranslation3d(-cyl(1:3));
0056 % Transform the points.
0057 ptTfm = transformPoint3d(ptProj, TFM);
0058 % Convert the transformed points to cylindrical coordinates.
0059 [ptsTheta, ptsRadius, ptsHeight] = cart2cyl(ptTfm);
0060 assert(ismembertol(ptsRadius(1),ptsRadius(2)))
0061 assert(ismembertol(ptsRadius(1),cylRadius))
0062 
0063 % Copy thetas for the conjugate geodesic
0064 ptsTheta(:,:,2) = ptsTheta;
0065 ptsTheta(1,1,2) = ptsTheta(1,1,2) + 2*pi;
0066 
0067 geoCyl = nan(n,3,size(ptsTheta,3));
0068 arcLength = nan(1,size(ptsTheta,3));
0069 for t = 1:size(ptsTheta,3)
0070     [geoCyl(:,:,t), arcLength(t)] = geoCurve(ptsTheta(:,:,t), cylRadius, ptsHeight, n);
0071 end
0072 
0073 % Select the shortest geodesic
0074 if arcLength(1) <= arcLength(2)
0075     % Transform the geodesics back to the global coordinate system
0076     geo = transformPoint3d(geoCyl(:,:,1), inv(TFM));
0077     conGeo = transformPoint3d(geoCyl(:,:,2), inv(TFM));
0078     geoLength = arcLength(1);
0079     conGeoLength = arcLength(2);
0080 else
0081     % Transform the geodesics back to the global coordinate system
0082     geo = transformPoint3d(geoCyl(:,:,2), inv(TFM));
0083     conGeo = transformPoint3d(geoCyl(:,:,1), inv(TFM));
0084     geoLength = arcLength(2);
0085     conGeoLength = arcLength(1);
0086 end
0087 
0088 end
0089 
0090 function [geo, arcLength] = geoCurve(theta, r, z, n)
0091 % Parametric expression of the geodesic curve
0092 u = linspace(theta(1),theta(2),n)';
0093 geo(:,1) = r*cos(u);
0094 geo(:,2) = r*sin(u);
0095 geo(:,3) = (z(2)-z(1))/(theta(2)-theta(1))*u + (z(1)*theta(2)-z(2)*theta(1))/(theta(2)-theta(1));
0096 if all(isnan(geo(:,3)))
0097     geo(:,3) = linspace(z(1),z(2),n)';
0098 end
0099 arcLength = sqrt(r^2*(theta(2)-theta(1))^2+(z(2)-z(1))^2);
0100 end