SurfStatEdg

Finds edges of a triangular mesh or a lattice.

 Usage: edg = SurfStatEdge( surf );

 surf.tri = t x 3 matrix of triangle indices, 1-based, t=#triangles.
 or
 surf.lat = 3D logical array, 1=in, 0=out.

 edg = e x 2 matrix of edge indices, 1-based, e=#edges.

0001 function edg = SurfStatEdge( surf )
0002 
0003 %Finds edges of a triangular mesh or a lattice.
0004 %
0005 % Usage: edg = SurfStatEdge( surf );
0006 %
0007 % surf.tri = t x 3 matrix of triangle indices, 1-based, t=#triangles.
0008 % or
0009 % surf.lat = 3D logical array, 1=in, 0=out.
0010 %
0011 % edg = e x 2 matrix of edge indices, 1-based, e=#edges.
0012 
0013 if isfield(surf,'tri')
0014     tri=sort(surf.tri,2);
0015     edg=unique([tri(:,[1 2]); tri(:,[1 3]); tri(:,[2 3])],'rows');
0016 end
0017 if isfield(surf,'lat')
0018     % See the comments of SurfStatResels for a full explanation:
0019     [I,J,K]=size(surf.lat);
0020     IJ=I*J;
0021     [i,j]=ndgrid(1:int32(I),1:int32(J));
0022     i=i(:);
0023     j=j(:);
0024     n1=(I-1)*(J-1)*6+(I-1)*3+(J-1)*3+1;
0025     n2=(I-1)*(J-1)*3+(I-1)+(J-1);
0026     edg=zeros((K-1)*n1+n2,2,'int32');
0027     for f=0:1
0028         c1=int32(find(rem(i+j,2)==f & i<I & j<J));
0029         c2=int32(find(rem(i+j,2)==f & i>1 & j<J));
0030         c11=int32(find(rem(i+j,2)==f & i==I & j<J));
0031         c21=int32(find(rem(i+j,2)==f & i==I & j>1));
0032         c12=int32(find(rem(i+j,2)==f & i<I & j==J));
0033         c22=int32(find(rem(i+j,2)==f & i>1 & j==J));
0034         edg0=[c1     c1+1; % bottom slice
0035             c1     c1+I;
0036             c1     c1+1+I;
0037             c2-1   c2;
0038             c2-1   c2-1+I;
0039             c2     c2-1+I;
0040             c11    c11+I;
0041             c21-I  c21;
0042             c12    c12+1;
0043             c22-1  c22];
0044         edg1=[c1     c1+IJ; % between slices
0045             c1     c1+1+IJ;
0046             c1     c1+I+IJ;
0047             c11    c11+IJ;
0048             c11    c11+I+IJ;
0049             c12    c12+IJ;
0050             c12    c12+1+IJ];
0051         edg2=[c2-1   c2-1+IJ;
0052             c2     c2-1+IJ;
0053             c2-1+I c2-1+IJ;
0054             c21-I  c21-I+IJ;
0055             c21    c21-I+IJ;
0056             c22-1  c22-1+IJ;
0057             c22    c22-1+IJ];
0058         if f
0059             for k=2:2:(K-1)
0060                 edg((k-1)*n1+(1:n1),:)=[edg0; edg2; edg1; IJ 2*IJ]+(k-1)*IJ;
0061             end
0062         else
0063             for k=1:2:(K-1)
0064                 edg((k-1)*n1+(1:n1),:)=[edg0; edg1; edg2; IJ 2*IJ]+(k-1)*IJ;
0065             end
0066         end
0067         if rem(K+1,2)==f
0068             edg((K-1)*n1+(1:n2),:)=edg0(1:n2,:)+(K-1)*IJ; % top slice
0069         end
0070     end
0071     % index by voxels in the lat
0072     vid=int32(cumsum(surf.lat(:)).*surf.lat(:));
0073     % only inside the lat
0074     edg=vid(edg(all(surf.lat(edg),2),:));
0075 end
0076 
0077 return
0078 end