surf2stl

SURF2STL   Write STL file from surface data.
   SURF2STL('filename',X,Y,Z) writes a stereolithography (STL) file
   for a surface with geometry defined by three matrix arguments, X, Y
   and Z.  X, Y and Z must be two-dimensional arrays with the same size.

   SURF2STL('filename',x,y,Z), uses two vector arguments replacing
   the first two matrix arguments, which must have length(x) = n and
   length(y) = m where [m,n] = size(Z).  Note that x corresponds to
   the columns of Z and y corresponds to the rows.

   SURF2STL('filename',dx,dy,Z) uses scalar values of dx and dy to
   specify the x and y spacing between grid points.

   SURF2STL(...,'mode') may be used to specify the output format.

     'binary' - writes in STL binary format (default)
     'ascii'  - writes in STL ASCII format

   Example:

     surf2stl('test.stl',1,1,peaks);

   See also SURF.

   Author: Bill McDonald, 02-20-04

0001 function surf2stl(filename,x,y,z,mode,writemode)
0002 %SURF2STL   Write STL file from surface data.
0003 %   SURF2STL('filename',X,Y,Z) writes a stereolithography (STL) file
0004 %   for a surface with geometry defined by three matrix arguments, X, Y
0005 %   and Z.  X, Y and Z must be two-dimensional arrays with the same size.
0006 %
0007 %   SURF2STL('filename',x,y,Z), uses two vector arguments replacing
0008 %   the first two matrix arguments, which must have length(x) = n and
0009 %   length(y) = m where [m,n] = size(Z).  Note that x corresponds to
0010 %   the columns of Z and y corresponds to the rows.
0011 %
0012 %   SURF2STL('filename',dx,dy,Z) uses scalar values of dx and dy to
0013 %   specify the x and y spacing between grid points.
0014 %
0015 %   SURF2STL(...,'mode') may be used to specify the output format.
0016 %
0017 %     'binary' - writes in STL binary format (default)
0018 %     'ascii'  - writes in STL ASCII format
0019 %
0020 %   Example:
0021 %
0022 %     surf2stl('test.stl',1,1,peaks);
0023 %
0024 %   See also SURF.
0025 %
0026 %   Author: Bill McDonald, 02-20-04
0027 
0028 error(nargchk(4,6,nargin));
0029 
0030 if (ischar(filename)==0)
0031     error( 'Invalid filename');
0032 end
0033 
0034 if (nargin < 5)
0035     mode = 'binary';
0036 elseif (strcmp(mode,'ascii')==0)
0037     mode = 'binary';
0038 end
0039 
0040 if (ndims(z) ~= 2)
0041     error( 'Variable z must be a 2-dimensional array' );
0042 end
0043 
0044 if any( (size(x)~=size(z)) | (size(y)~=size(z)) )
0045     
0046     % size of x or y does not match size of z
0047     
0048     if ( (length(x)==1) & (length(y)==1) )
0049         % Must be specifying dx and dy, so make vectors
0050         dx = x;
0051         dy = y;
0052         x = ((1:size(z,2))-1)*dx;
0053         y = ((1:size(z,1))-1)*dy;
0054     end
0055         
0056     if ( (length(x)==size(z,2)) & (length(y)==size(z,1)) )
0057         % Must be specifying vectors
0058         xvec=x;
0059         yvec=y;
0060         [x,y]=meshgrid(xvec,yvec);
0061     else
0062         error('Unable to resolve x and y variables');
0063     end
0064         
0065 end
0066 
0067 if strcmp(mode,'ascii')
0068     % Open for writing in ascii mode
0069    % fid = fopen(filename,'w');
0070    fid = fopen(filename,writemode);
0071 else
0072     % Open for writing in binary mode
0073     fid = fopen(filename,'wb+');
0074 end
0075 
0076 if (fid == -1)
0077     error( sprintf('Unable to write to %s',filename) );
0078 end
0079 
0080 title_str = sprintf('Created by surf2stl.m %s',datestr(now));
0081 
0082 if strcmp(mode,'ascii')
0083     fprintf(fid,'solid %s\r\n',title_str);
0084 else
0085     str = sprintf('%-80s',title_str);    
0086     fwrite(fid,str,'uchar');         % Title
0087     fwrite(fid,0,'int32');           % Number of facets, zero for now
0088 end
0089 
0090 nfacets = 0;
0091 
0092 for i=1:(size(z,1)-1)
0093     for j=1:(size(z,2)-1)
0094         
0095         p1 = [x(i,j)     y(i,j)     z(i,j)];
0096         p2 = [x(i,j+1)   y(i,j+1)   z(i,j+1)];
0097         p3 = [x(i+1,j+1) y(i+1,j+1) z(i+1,j+1)];
0098         val = local_write_facet(fid,p1,p2,p3,mode);
0099         nfacets = nfacets + val;
0100         
0101         p1 = [x(i+1,j+1) y(i+1,j+1) z(i+1,j+1)];
0102         p2 = [x(i+1,j)   y(i+1,j)   z(i+1,j)];
0103         p3 = [x(i,j)     y(i,j)     z(i,j)];        
0104         val = local_write_facet(fid,p1,p2,p3,mode);
0105         nfacets = nfacets + val;
0106         
0107     end
0108 end
0109 
0110 if strcmp(mode,'ascii')
0111     fprintf(fid,'endsolid %s\r\n',title_str);
0112 else
0113     fseek(fid,0,'bof');
0114     fseek(fid,80,'bof');
0115     fwrite(fid,nfacets,'int32');
0116 end
0117 
0118 fclose(fid);
0119 
0120 disp( sprintf('Wrote %d facets',nfacets) );
0121 
0122 
0123 
0124 % Local subfunctions
0125 
0126 function num = local_write_facet(fid,p1,p2,p3,mode)
0127 
0128 if any( isnan(p1) | isnan(p2) | isnan(p3) )
0129     num = 0;
0130     return;
0131 else
0132     num = 1;
0133     n = local_find_normal(p1,p2,p3);
0134     
0135     if strcmp(mode,'ascii')
0136         
0137         fprintf(fid,'facet normal %.7E %.7E %.7E\r\n', n(1),n(2),n(3) );
0138         fprintf(fid,'outer loop\r\n');        
0139         fprintf(fid,'vertex %.7E %.7E %.7E\r\n', p1);
0140         fprintf(fid,'vertex %.7E %.7E %.7E\r\n', p2);
0141         fprintf(fid,'vertex %.7E %.7E %.7E\r\n', p3);
0142         fprintf(fid,'endloop\r\n');
0143         fprintf(fid,'endfacet\r\n');
0144         
0145     else
0146         
0147         fwrite(fid,n,'float32');
0148         fwrite(fid,p1,'float32');
0149         fwrite(fid,p2,'float32');
0150         fwrite(fid,p3,'float32');
0151         fwrite(fid,0,'int16');  % unused
0152         
0153     end
0154     
0155 end
0156 
0157 
0158 function n = local_find_normal(p1,p2,p3)
0159 
0160 v1 = p2-p1;
0161 v2 = p3-p1;
0162 v3 = cross(v1,v2);
0163 n = v3 ./ sqrt(sum(v3.*v3));