make_cylinder_plot

 Make cylinder around set of points which represent sequential line
  [V,F,X] = make_cylinder(X,dN,dr)
 X  : set of points 3D coordinates [N x 3]
 dN : number of points consist cylinder
      = 6 means hexagonal cylinder will be made
 dr : radius of cylinder (same unit with X)
 V  : vertex coordinate for cylindrical lines
 F  : patch index for cylindrical lines

 2014-11-10 Masa-aki Sato

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function    [V,F,X,Info] = make_cylinder(X,dN,dr)
0002 % Make cylinder around set of points which represent sequential line
0003 %  [V,F,X] = make_cylinder(X,dN,dr)
0004 % X  : set of points 3D coordinates [N x 3]
0005 % dN : number of points consist cylinder
0006 %      = 6 means hexagonal cylinder will be made
0007 % dr : radius of cylinder (same unit with X)
0008 % V  : vertex coordinate for cylindrical lines
0009 % F  : patch index for cylindrical lines
0010 %
0011 % 2014-11-10 Masa-aki Sato
0012 %
0013 % Copyright (C) 2011, ATR All Rights Reserved.
0014 % License : New BSD License(see VBMEG_LICENSE.txt)
0015 
0016 %　各線分に沿って dN 角柱を作り、patchを生成する
0017 
0018 [N, d] = size(X);
0019 if d ~= 3, error('data dimension error'); end
0020 if N < 2, error('There is only one point'); end
0021 
0022 % dx : line segment direction [N x 3]
0023 dx = diff(X);
0024 dx = [dx; dx(end,:)];
0025 
0026 % Find reverse direction
0027 dxs = sum(dx(1:N-1,:).*dx(2:N,:),2);
0028 dxs = [1; dxs];
0029 irv = find(dxs < 0);
0030 Nrev = length(irv);
0031 
0032 % Insert duplicate points at reverse points
0033 Xnew = [];
0034 ist = 1;
0035 for m = 1:Nrev
0036     n = irv(m);
0037     Xnew = [Xnew; X(ist:n,:); X(n,:)];
0038     ist = n+1;
0039 end
0040 
0041 if ist <= N
0042     Xnew = [Xnew; X(ist:N,:)];
0043 end
0044 
0045 if size(Xnew,1) ~= (N+Nrev), error('Insert error'); end;
0046 
0047 X = Xnew;
0048 N = size(X,1);
0049 
0050 % Recalculate line segment direction [N x 3]
0051 dx = diff(X);
0052 dx = [dx; dx(end,:)];
0053 
0054 % Reverse direction check
0055 dxs = sum(dx(1:N-1,:).*dx(2:N,:),2);
0056 dxs = [1; dxs];
0057 irv = find(dxs < 0);
0058 Nrev = length(irv);
0059 
0060 [dY,dZ] = make_orthogonal_vec(dx);
0061 
0062 % Find zero difference vector
0063 dxx = sum(dx.^2,2);
0064 zx  = find(dxx == 0);
0065 Nzero = length(zx);
0066 
0067 Info.Nzero = Nzero;
0068 Info.Nrev  = Nrev;
0069 
0070 % X(n-1) -> X(n) = X(n+1) -> X(n+2)
0071 %    dx(n-1)  dx(n)=0  dx(n+1)
0072 % X(n-1) -> X(n)  : stop line
0073 %             dx(n)=dx(n-1)
0074 %      start line : X(n+1) -> X(n+2)
0075 
0076 % Copy preceding direction to stop connected cylinder
0077 for m=1:Nzero
0078     n = zx(m);
0079     if n > 1
0080         dx(n,:) = dx(n-1,:);
0081         dY(n,:) = dY(n-1,:);
0082         dZ(n,:) = dZ(n-1,:);
0083     end
0084 end
0085 
0086 % incremental angle of cylindrical points
0087 ang = 2*pi/dN;
0088 da  = ang*(0:dN-1);
0089 % [Y,Z] coordinate of cylindrical points
0090 y = cos(da);
0091 z = sin(da);
0092 
0093 % make dN points for each N points
0094 V   = zeros(dN*N,3);
0095 xyz = zeros(N,3,dN);
0096 
0097 % make cylindrical points around X with radius dr
0098 for m=1:dN
0099     xyz(:,:,m) = X + (dY*y(m) + dZ*z(m))*dr;
0100 end
0101 
0102 % find nearest point from xyz(n,:,:) to xyz(n+1,:,:)
0103 % which is cylindrical points in next X(n+1,:)
0104 % and reorder the next cylindrical points xyz(n+1,:,:)
0105 for n=1:N-1
0106     % search for min distance between xyz(n,:,:) and xyz(n+1,:,:)
0107     xp = reshape(xyz(n,:,:),[3 dN]);
0108     xn = reshape(xyz(n+1,:,:),[3 dN]);
0109     dd = sq_distance(xp,xn);
0110     
0111     [dmin, jmin] = min(dd,[],2);
0112     [dtmp, kmin] = min(dmin);
0113     % min distance : xyz(n,:,kmin) - xyz(n+1,:,jmin(kmin))
0114     
0115     imin = jmin(kmin) - kmin + 1;
0116     
0117     % reorder the next cylindrical points xyz(n+1,:,:)
0118     ind = imin + (0:dN-1);
0119     ind = mod(ind-1,dN) + 1;
0120     
0121     xyz(n+1,:,:) = xyz(n+1,:,ind);
0122 end
0123 
0124 % cylindrical points of X(n,:)
0125 % : xyz(n,:,m)
0126 % = V( dN*(n-1) + m, :)
0127 nj = 0:dN:(dN*(N-1));
0128 
0129 for m=1:dN
0130     V(nj + m,:) = xyz(:,:,m);
0131 end
0132 
0133 %--- Example of vertex & patch
0134 %  1  2  3  1
0135 %  |/ |/ |/ |
0136 %  4  5  6  4
0137 %---  n=1, dN=3
0138 %  F = [n n+1 n+dN;
0139 %         n+1 n+dN n+1+dN]
0140 % There are 2*dN patches for pair of consecutive cylindrical points
0141 F  = zeros(2*dN*(N-1),3);
0142 
0143 % make patch index
0144 for i=1:N-1
0145     np = dN*(i-1);  % start point index
0146     nf = 2*dN*(i-1);% start patch index
0147     
0148     for n=1:dN
0149         if n==dN, 
0150             m = 1;
0151         else
0152             m = n+1;
0153         end
0154         
0155         %  F = [n n+1 n+dN;
0156         %         n+1 n+dN n+1+dN]
0157         F(nf+n,:)    = [n   m    n+dN] + np;
0158         F(nf+dN+n,:) = [m   m+dN n+dN] + np;
0159     end
0160 end