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);
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));
if max(max(abs(A-A'))) == 0
% found it
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;
% Assert that edges correspond to face sharers
% 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!