Consistent (tet) mesh neighbor information

Alec Jacobson

January 10, 2013

weblog/

I've been using tet mesh from a matlab interface I wrote. It produces a list of tet indices, I call TT in matlab, where each row is a four-tuple of tet indices. It can also output a .neigh file which contains a list of each tet's neighbor by index. I call this TN in matlab. This is also a list of rows of four, when now a positive index of j on row i means that tet i is a neighbor with tet j. Neighboring means sharing a face. Sharing a face means tet i and tet j have 3 indices in common. One would expect that these come in a consistent ordering. For example, if j appears as the kth element on row i, then I might expect tet j to be the "kth neighbor" of tet i. So if k=3 I expect that tet i shares the face opposite its kth vertex (kth column in TT). ~~Doch leider ist es nicht wahr.~~ (actually it is, see below). Instead it seemed that tetgen was using a different ordering. To find out which consistent ordering scheme it was using I wrote a little script:

TE = [repmat(1:size(TT,1),1,size(TT,2))' TN(:)];
% verify edges are symmetric
A = sparse(TE(TE(:,2)>0,1),TE(TE(:,2)>0,2),1);
assert(size(A,1)==size(A,2));
assert(max(max(A'-A))==0);

all_perms = perms(1:4);
% loop over all perms
for p = 1:size(all_perms,1)
  TI = [T(:,[2 3 4]); T(:,[3 4 1]); T(:,[4 1 2]); T(:,[1 2 3])];
  TI = [ ...
    T(:,all_perms(p,[2 3 4])); ...
    T(:,all_perms(p,[3 4 1])); ...
    T(:,all_perms(p,[4 1 2])); ...
    T(:,all_perms(p,[1 2 3]))];
  A = sparse(TE(TE(:,2)>0,1),TE(TE(:,2)>0,2),sum(TI(TE(:,2)>0),2));
  fprintf('%d: %d\n',p,full(max(max(abs(A-A')))));
  if max(max(abs(A-A'))) == 0
    % found it
    break;
  end
end

First I do a little sanity check to make sure that tet's agree that they are neighbors. Then I consider every (4! = 24) possible consistent ordering of tets by checking whether their shared faces agree. But again! Foiled. Tetgen seems to output neighbors in an inconsistent manner. To resolve this I wrote another small script:

% Get rid of boundary edges
TE = TE(TE(:,2)>0,:);
TNI = zeros(size(TE,1),1);
for c = 1:size(TT,2)
  [TNIc] = all(~bsxfun(@eq,TT(TE(:,1),c),TT(TE(:,2),:)),2);
  % Assert face-manifoldness
  assert(all(TNI(TNIc) == 0));
  TNI(TNIc) = c;
end
% Assert that edges correspond to face sharers
assert(all(TNI>0));
% Rebuild TN
new_TN = full(sparse(TE(:,1),TNI,TE(:,2),size(TN,1),size(TN,2)));
new_TN(new_TN==0) = -1;

As a consequence the first script can be used to check for a consistent ordering and the second script can be used to check for face-manifoldness and correct neighbor relations. Update: Turns out tetgen was fine and I had a bug. I was reordering the columns in my TT, without reordering the columns in TN. But these scripts did help me find and eliminate my bug!