extrudecurve

 [node,face,yz0,yz1]=extrudecurve(xy, yz, Nx, Nz, Nextrap, spacing, anchor)
 
 create a triangular surface mesh by swining a 2D spline along another 2D
 spline curve

 author: Qianqian Fang, <q.fang at neu.edu>

 input:
      xy: a 2D spline path, along which the surface is extruded, defined
          on the x-y plane
      yz: a 2D spline which will move along the path to form a surface,
          defined on the y-z plane
      Nx: the count of sample points along the extrusion path (xy), if
          ignored, it is 40
      Nz: the count of sample points along the curve to be extruded (yz),
          if ignored, it is 40
      Nextrap: number of points to extrapolate outside of the xy/yz
          curves, 0 if ignored
      spacing: define a spacing scaling factor for spline interpolations,
          1 if ignored
      anchor: the 3D point in the extruded curve plane (yz) that is aligned
          at the nodes long the extrusion path. this point does not have
          to be located on the yz curve; orig = [ox oy oz], if ignored, it
          is set as the point on the interpolated yz with the largested
          y-value
      dotopbottom: a flag, if set to 1, tessellated top and bottom faces
          will be added. default is 0.

 output:
      node: 3D node coordinates for the generated surface mesh
      face: triangular face patches of the generated surface mesh, each
           row represents a triangle denoted by the indices of the 3 nodes

 -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)

0001 function [node,face,yz0,yz1]=extrudecurve(xy, yz, Nx, Nz, Nextrap, spacing, anchor, dotopbottom)
0002 %
0003 % [node,face,yz0,yz1]=extrudecurve(xy, yz, Nx, Nz, Nextrap, spacing, anchor)
0004 %
0005 % create a triangular surface mesh by swining a 2D spline along another 2D
0006 % spline curve
0007 %
0008 % author: Qianqian Fang, <q.fang at neu.edu>
0009 %
0010 % input:
0011 %      xy: a 2D spline path, along which the surface is extruded, defined
0012 %          on the x-y plane
0013 %      yz: a 2D spline which will move along the path to form a surface,
0014 %          defined on the y-z plane
0015 %      Nx: the count of sample points along the extrusion path (xy), if
0016 %          ignored, it is 40
0017 %      Nz: the count of sample points along the curve to be extruded (yz),
0018 %          if ignored, it is 40
0019 %      Nextrap: number of points to extrapolate outside of the xy/yz
0020 %          curves, 0 if ignored
0021 %      spacing: define a spacing scaling factor for spline interpolations,
0022 %          1 if ignored
0023 %      anchor: the 3D point in the extruded curve plane (yz) that is aligned
0024 %          at the nodes long the extrusion path. this point does not have
0025 %          to be located on the yz curve; orig = [ox oy oz], if ignored, it
0026 %          is set as the point on the interpolated yz with the largested
0027 %          y-value
0028 %      dotopbottom: a flag, if set to 1, tessellated top and bottom faces
0029 %          will be added. default is 0.
0030 %
0031 % output:
0032 %      node: 3D node coordinates for the generated surface mesh
0033 %      face: triangular face patches of the generated surface mesh, each
0034 %           row represents a triangle denoted by the indices of the 3 nodes
0035 %
0036 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0037 %
0038 
0039 if(nargin<3)
0040     Nx=30;
0041 end
0042 if(nargin<4)
0043     Nz=30;
0044 end
0045 if(nargin<5)
0046     Nextrap=0;
0047 end
0048 if(nargin<6)
0049     spacing=1;
0050 end
0051 
0052 xrange=max(xy(:,1))-min(xy(:,1));
0053 dx=xrange/Nx;
0054 xi=min(xy(:,1))-Nextrap*dx:spacing*dx:max(xy(:,1))+Nextrap*dx;
0055 pxy = spline(xy(:,1), xy(:,2));
0056 
0057 yi = ppval(pxy,xi);
0058 dy = gradient(yi);
0059 dx = gradient(xi);
0060 
0061 nn=sqrt(dx.*dx + dy.*dy);
0062 normaldir=[dx(:)./nn(:) dy(:)./nn(:)];
0063 
0064 zrange=max(yz(:,2))-min(yz(:,2));
0065 dz=zrange/Nz;
0066 zi=min(yz(:,2))-Nextrap*dz:spacing*dz:max(yz(:,2))+Nextrap*dz;
0067 pyz = spline(yz(:,2), yz(:,1));
0068 
0069 yyi = ppval(pyz,zi);
0070 
0071 if(~exist('anchor','var') || isempty(anchor))
0072     [ymax, loc]=max(yyi);
0073     anchor=[0 yyi(loc) zi(loc)];
0074 end
0075     
0076 node=zeros(length(zi)*length(xi),3);
0077 face=zeros(2*(length(zi)-1)*(length(xi)-1),3);
0078 
0079 xyz=[yyi(:) yyi(:) zi(:)];
0080 xyz(:,1)=0;
0081 for i=1:length(xi)
0082     rot2d=[normaldir(i,1) -normaldir(i,2); normaldir(i,2) normaldir(i,1)];
0083     offset=[xi(i) yi(i) anchor(2)];
0084     newyz=[((rot2d*(xyz(:,1:2)-repmat(anchor(1:2),size(xyz,1),1))' + repmat(offset(:,1:2),length(zi),1)'))' , xyz(:,end)];
0085     node((i-1)*length(zi)+1:i*length(zi),:)=newyz;
0086     if(i>1)
0087        face((i-2)*2*(length(zi)-1)+1:(i-1)*2*(length(zi)-1),:)= ...
0088                        [ [1:length(zi)-1 ;  (1-length(zi):-1) ; 2:length(zi)]+(i-1)*length(zi)  ...
0089                        [ (1-length(zi):-1); (2-length(zi):0)  ; 2:length(zi)]+(i-1)*length(zi) ]';
0090     end
0091     if(i==Nextrap+1)
0092         yz0=newyz(Nextrap+1:end-Nextrap,:);
0093     end
0094     if(i==length(xi)-Nextrap)
0095         yz1=newyz(Nextrap+1:end-Nextrap,:);
0096     end
0097 end
0098 
0099 % add two flat polygons on the top and bottom of the contours
0100 % to ensure the enclosed surface is not truncated by meshfix
0101 
0102 if(nargin>=8 && dotopbottom==1)
0103     nump = length(xi);
0104     C = [(1:(nump-1))' (2:nump)'; nump 1];
0105     dt = delaunayTriangulation(xi(:), yi(:), C);
0106     io = dt.isInterior();
0107     endface=dt(io,:);
0108     endface=(endface-1)*size(newyz,1)+1;
0109 
0110     % append the top/bottom faces to the extruded mesh
0111     face=[face; endface; endface+size(newyz,1)-1];
0112 end
0113 
0114 [node,face]=meshcheckrepair(node,face,'dup');