function b = gmoffread(filename, tol)% gmoffread: read a brep from a file in OFF format% b = gmoffread(filename)% This routine reads a brep in OFF format from a file.% See the documentation for a description of OFF format.% OFF format was developed at the University of Minnesota% Geometry Center.if nargin < 2  global GM_DEFAULT_TOL;  tol = GM_DEFAULT_TOL;endfi = fopen(filename,'r');if fi < 0  error(sprintf('File %s not found', filename))endlineco = 1;line = gm_strim(fgetl(fi));lineco = lineco + 1;if ~strcmpi(line, 'off')  error('File must begin with keyword OFF')enda = sscanf(fgetl(fi), '%d %d %d');lineco = lineco + 1;if length(a) ~= 3  error ('Second line of file must contain three integers')endnumvtx = a(1);numfacet = a(2);numedge = a(3);cplist = zeros(3,numvtx);topvertices = cell(5,numvtx);for vnum = 1 : numvtx  coord = sscanf(fgetl(fi), '%e %e %e');  if length(coord) ~= 3    errmsg = sprintf('Each vertex line must contain three coordinates (error on line %d)', lineco);    error(errmsg)  end  lineco = lineco + 1;  cplist(:,vnum) = coord;  vlist{1,vnum} = sprintf('v%d', vnum - 1);  vlist{2,vnum} = {};  vlist{3,vnum} = {};  vlist{4,vnum} = {};  vlist{5,vnum} = {'vertex'; []; [vnum - 1]};endelist = {};slist = cell(5,numfacet);basicscanstr = '%e ';ecount = sparse([],[],[],numvtx,numvtx);ehash = sparse([],[],[],numvtx,numvtx);totalnume = 0;rbdry = {};for fnum = 1 : numfacet  li = fgets(fi);  numv = sscanf(li, '%d');  if length(numv) == 0    errmsg = sprintf('Missing vertex count field in facet entry (line #%d)', ...                       lineco);    error(errmsg);  end  numv = numv(1);  idx = [1;2;3] * ones(1,numv + 4);  scanstring = basicscanstr(idx(:));  rval = sscanf(li, scanstring);  if length(rval) ~= numv + 1 & length(rval) ~= numv + 4    errmsg = sprintf('Facet has wrong number of entries (line #%d)', lineco);    error(errmsg);  end  if length(rval) == numv + 1    color = [];  else    color = rval(numv + 2 : numv + 4);    if any(color < 0) | any(color > 1)      error(sprintf('Color entry out of range (must be in [0,1]) (line #%d)', ...          lineco));    end  end  rval = rval(2:numv+1);% deduplicate rval  while any(rval(1:numv-1) == rval(2:numv))    f1 = find(rval(1:numv-1) == rval(2:numv));    f1 = f1(1);    rval = [rval(1:f1);rval(f1+2:end)];    numv = numv - 1;  end     if any(rval < 0 | rval >= numvtx)    error(sprintf('Vertex index out of range (must be in 0..%d) (line #%d)', ...          numvtx - 1, lineco));  end  lineco = lineco + 1;  sbdry = cell(1,numv);  thisvlist = zba(zeros(numv,3));  thiselist = zeros(numv,2);  for k = 1 : numv    v1 = rval(k);    if k == numv      nextk = 1;    else      nextk = k + 1;    end    v2 = rval(nextk);    thisvlist(k-1,:) = cplist(:,v1 + 1)';    thiselist(k,:) = [k - 1,nextk - 1];    if v2 < v1      tmp = v2;      v2 = v1;      v1 = tmp;    end    ecount(v1+1,v2+1) = ecount(v1+1,v2+1) + 1;    topidx = ehash(v1 + 1, v2 + 1);    if topidx == 0      totalnume = totalnume + 1;      elist{1,totalnume} = sprintf('e%d', totalnume - 1);      elist{2,totalnume} = {};      elist{3,totalnume} = {sprintf('v%d',v1),sprintf('v%d',v2)};      elist{4,totalnume} = {};      elist{5,totalnume} = {'bezier_curve';[1];[v1,v2]};      ehash(v1+1,v2+1) = totalnume;      topidx = totalnume;    end    sbdry{k} = sprintf('e%d', topidx - 1);  end  trilist = gm_polytri(thisvlist, thiselist, tol);  [numtri,scrap] = size(trilist);  bezlist = cell(3,numtri);  for l = 1 : numtri    bezlist{1,l} = 'bezier_triangle';    bezlist{2,l} = [1];    bezlist{3,l} = rval(trilist(l,:)+[1,1,1]);  end   slist{1,fnum} = sprintf('s%d', fnum - 1);  rbdry{length(rbdry)+1} = sprintf('s%d', fnum - 1);  if length(color) == 0    slist{2,fnum} = {};  else    slist{2,fnum} = {'color'; sprintf('(%f %f %f 1)', ...      color(1),color(2),color(3))};  end  slist{3,fnum} = sbdry;  slist{4,fnum} = {};  slist{5,fnum} = bezlist;end[el1,el2, v] = find(ecount);emit2d = 0;incorrecte = find(v~=2);for jj = 1 : length(incorrecte)  x = incorrecte(jj);  disp(sprintf('Boundary does not close: edge %d,%d is contained by %d polygon(s) (should be 2)', el1(x), el2(x), v(x)))  emit2d = 1;endglobal GM_BREP_TYPE_CODEif emit2d == 1  disp('Because boundary does not close, returned brep will have no regions');  b = zba({GM_BREP_TYPE_CODE; 2; 3; {}; cplist; vlist; elist; slist});else  rlist = {'region'; {}; rbdry; {}; {}};  b = zba({GM_BREP_TYPE_CODE; 3; 3; {}; cplist; vlist; elist; slist; rlist});end% ------------------------------------------------------------------% Copyright (c) 1999 by Cornell University.  All rights reserved.% See the accompanying file 'Copyright' for authorship information,% the terms of the license governing this software, and disclaimers% concerning this software.% ------------------------------------------------------------------% This file is part of the QMG software.  % Version 2.0 of QMG, release date September 3, 1999% ------------------------------------------------------------------