ctri_to_mesh_converter.py

import numpy as np

import sys
from collections import defaultdict


"""
Converter for CGAL's  tetrahedral file format.

This requires a bunch of routines to ensure that the tetrahedra
are correctly oriented.
This includes reorienting "infinite" tetradedra that have an 
OOB vertex / vertx at infinity.
The orientation of these infinite tets is defined as the opposite
of finite mirror.

What is a "mirror"? Every finite 2D face (triangle) is shared by exactly 2 3D tets. These tets are "mirrors" of each other allong this 2D face.
In 2D (with 1D mirrors), the picture looks like this:
a----b
|   /|
|1 / |
| /  |
|/  2|
c----d
Triangle abc and dbc share edge bc.
Triangle acb is counter-clockwise (CCW); 
its mirrored version dcb is CW.
All tets must be in CCW orientation (e.g. cdb).

If the other file format does not specify these infinite tets, but
does specify boundary faces (like INRIA medit files), then these
infinite tets can be generated by adding the oob vertex to the node list
and swapping to ensure correct orientation.

The file also requires adjacency information for every tet.
So if a tet is abcd, then we need to generate a list of 4 indices
refering to mirror tetrahedra for every face.
The order is defined by dropping vertex indices in the propertly oriented tet.
E.g. mirrors along the faces [bcd,acd,abd,abc].
"""

##############################################
# File reading stuff

# Strip comments, clean, and flip (so we can pop lines)
def prep_input(lines):
    lines = [line.split('#',1)[0].strip() for line in lines] # Strip comments
    lines = [line for line in lines if line] # Remove empty
    lines = list(reversed(lines)) # reverse for popping
    return lines

"""
Read in the output from CGAL::Triangulation_3's out stream operator (<<)
"""
def read_cgal_tri(filename):
    fh = open(filename,'r')
    lines = fh.readlines()
    fh.close()
    lines = prep_input(lines)

    vertices    = []
    tetrahedra = []
    neighbors = []

    D = int(lines.pop())
    V = int(lines.pop())
    vertices.append([np.nan]*3) # Explicitly add the oob node
    for v in xrange(V):
        coords = map(float,lines.pop().split())
        assert D == len(coords)
        vertices.append(coords)
    T = int(lines.pop())
    for t in xrange(T):
        idx = map(int,lines.pop().split())
        assert (D+1) == len(idx)
        tetrahedra.append(list(sorted(idx)))
    a = 0
    while lines:
        idx = map(int,lines.pop().split())
        assert (D+1) == len(idx)
        neighbors.append(list(sorted(idx)))
        a += 1
    return vertices,tetrahedra,neighbors

def read_inria_mesh(fileName):
    FH = open(fileName,'r')
    lines = FH.readlines()
    FH.close()
    lines = prep_input(lines)

    I = 0
    names = ['vertices','triangles','tetrahedra']
    types = [float,int,int]
    objects = {}
    for (name,T) in zip(names,types):
        lines_copy = list(lines)
        line = lines_copy.pop()
        while name not in line.lower():
            line = lines_copy.pop()
        assert name in line.lower()
        n = int(lines_copy.pop())
        # Ignore boundary information
        objs = []
        for i in xrange(n):
            line = lines_copy.pop()
            obj = map(T, line.split()[:-1])
            objs.append(obj)
        exec(name + ' = objs')
    vertices = [[np.nan]*3] + vertices # prepend oob
    triangles = [list(sorted(tri)) for tri in triangles]
    tetrahedra = [list(sorted(tet)) for tet in tetrahedra]
     
    return vertices,triangles,tetrahedra

########################################33
# File writing stuff

"""
Write information to a file that can be read into CGAL::Triangulation_3 using
the stream operator (>>)
NB: The tetrahedra MUST be in CCW order---see the reorienting function.
In particular, any non-finite tet's (tetrahedra that contain the 0 vertex id)
must have a mirror tet in CW order; see the mirror tet routines for 
more information
"""
def write_cgal_tri(filename,vertices,tetrahedra,neighbors):
    print "Writing to CGAL::TRI file"
    print "\tVertices",len(vertices)
    print "\tTetrahedra",len(tetrahedra)
    V = len(vertices)
    assert V > 0
    D = len(vertices[0])
    T = len(tetrahedra)
    
    SS = []
    SS.append(D)
    SS.append(V-1)
    SS.append(vertices[1:]) # Don't write OOB node
    SS.append(T)
    SS.append(tetrahedra)
    SS.append(neighbors)

    S = []
    for x in SS:
        if isinstance(x,list):
            for y in x:
                S.append(' '.join(map(str,y)))
        else:
            S.append(str(x))

    fh = open(filename,'w')    
    fh.write('\n'.join(S))
    fh.close()

def write_inria_mesh(filename,vertices,triangles,tetrahedra):
    print "Writing to INRIA mesh"
    print "\tVertices",len(vertices)
    print "\tTetrahedra",len(tetrahedra)
    names = ['Vertices','Triangles','Tetrahedra']
    SS = []
    SS.append("MeshVersionFormatted 1") # ?
    SS.extend(["Dimension",3])
    
    vertices = vertices[1:] # Remove oob
    
    for (name,objs) in zip(names,[vertices,triangles,tetrahedra]):
        SS.extend([name,len(objs)])
        for obj in objs:
            SS.append('\t'.join(map(str,obj + [1])))
            # Add uninformative boundary marker
    SS.append('End')

    SS = map(str,SS)
    
    fh = open(filename,'w')    
    fh.write('\n'.join(SS))
    fh.close()  

##############################################3
# Tetrahedron orientation

"""
Get face hashs from a tet
"""
face_hash = lambda x,i: tuple(sorted(x[:i] + x[(i+1):]))
def face_keys(tet):
    assert 4 == len(tet)
    keys = []
    for i in xrange(4):
        key = face_hash(tet,i)
        keys.append(key)
    return keys
    
"""
Turns 4 tetrahedra indices into a 4x3 matrix of
physical coordinates
"""
def tetrahedra_vertex_matrix(tet_idx,vertices):
    tet_list = [vertices[i] for i in tet_idx]
    tet_mat = np.array(tet_list)
    assert (4,3) == tet_mat.shape
    return tet_mat

def orient3d(tet_mat):
    assert (4,3) == tet_mat.shape
    assert not np.any(np.isnan(tet_mat))

    # Get difference vectors from the first vertex
    diff = tet_mat[1:,:] - tet_mat[0,:]
    assert (3,3) == diff.shape

    # The sign of the difference matrix is
    # the orientation
    orient = np.linalg.det(diff)
    return np.sign(orient)

def reorient_tetrahedras(vertices,tetrahedra,mirror):
    for (tet_id,tet) in enumerate(tetrahedra):
        if 0 in tet:
            # Mirror world            
            mirror_tet = get_mirror_tetrahedra(mirror,tetrahedra,tet_id,0)
            tet_mat = tetrahedra_vertex_matrix(mirror_tet,vertices)
            sign = -1
        else:
            # Regular world.
            sign = 1
            tet_mat = tetrahedra_vertex_matrix(tet,vertices)
        assert (4,3) == tet_mat.shape
        assert not np.any(np.isnan(tet_mat))
        orientation = orient3d(tet_mat)
        if orientation == 0:
            print "Degenerate tetrahedron:",tet_id,tet
            assert orientation != 0
        if sign*orientation < 0:
            # Flip. All even permutations are equivalent
            tetrahedra[tet_id][0],tetrahedra[tet_id][1] \
                = tetrahedra[tet_id][1],tetrahedra[tet_id][0]

####################################################3
# Functions for mirror and neighbor structure
            
"""
Every finite face (triangle) belongs to two tetrehedra.
The mirror data structure keeps track of this.
"""
def build_mirror(tetrahedra):
    mirror = defaultdict(set)
    for (tet_id,tet) in enumerate(tetrahedra):
        assert 4 == len(tet)
        keys = face_keys(tet)
        for key in keys:
            mirror[key].add(tet_id)
            if len(mirror[key]) > 2:
                print 'Triangle', key,\
                    'associated with too many tetrahedra:',\
                    mirror[key]
                assert len(mirror[key]) <= 2
    return mirror

"""
Takes a tet and a vertex id; 
Returns the mirror tet (vertex id list) for the vertex that is different
The order of the mirror tet list is the SAME as the original list,
but the different vertices have been swapped.
Therefore their orientation will be opposite for any finite tet pair.
For an infinite-finite tet pair, the infinite tet's orientation will be
DEFINED as the opposite of the finite tet's.
"""
def get_mirror_tetrahedra(mirror,tetrahedra,tet_id,vert_id):
    # Get the vertex id list
    tet = list(tetrahedra[tet_id])
    assert 4 == len(tet)
    assert vert_id in tet
    
    # Find the index, in the tet list of the supplied vertex id
    idx = tet.index(vert_id)
    
    # Get the face of the tet that doesn't use this vertex
    # This is the 'mirror'
    key = face_hash(tet,idx)

    # Get the pair of tet indices that share that share the mirror
    twins = mirror[key]
    assert tet_id in twins

    # Get the mirror vertex that is not in the face
    mirror_tet_id =  twins - set([tet_id])
    assert 1 == len(mirror_tet_id)
    mirror_tet_id = list(mirror_tet_id)[0] # Get unique tet index
    mirror_tet = tetrahedra[mirror_tet_id] # Get the vertex id list

    # What is different between the original vertex id list
    # and the mirror list?
    mirror_vert_id = set(mirror_tet) - set(tet)
    assert 1 == len(mirror_vert_id) 
    mirror_vert_id = list(mirror_vert_id)[0] # Get unique difference

    # Mirror tet vertex id list,
    # In same order as original tet vertex list, but with
    # the indicated vert replaced.
    # Will have the opposite orientation as the original.
    mirror_tet_replace = tet
    mirror_tet_replace[idx] = mirror_vert_id

    # Same vertex ids; perhaps different order
    assert(list(sorted(mirror_tet_replace)) == list(sorted(mirror_tet)))
    
    return mirror_tet_replace

"""
Generate the neighbors for each tetrahedra
Needs to be done in a specific order.
If the tet is abcd, then the first neighbor must be the
mirror tet that shares face bcd (e.g. is any even perturbation of b*cd)
"""
def generate_neighbors(tetrahedra,mirror):
    neighbors = []
    for (tet_id,tet) in enumerate(tetrahedra):
        tet_adj = []
        keys = face_keys(tet)
        for (face_id,key) in enumerate(keys):
            assert tet_id in mirror[key]
            mirror_tet_id = mirror[key] - set([tet_id])
            assert 1 == len(mirror_tet_id)
            mirror_tet_id = list(mirror_tet_id)[0]
            tet_adj.append(mirror_tet_id)
        neighbors.append(tet_adj)
    return neighbors


"""
Do triangles have a "correct orientation" in 3D?
Is there a "righthand rule" here?
"""
def generate_boundary_triangles(vertices,tetrahedra):
    triangles = []
    for tet in tetrahedra:
        if 0 not in tet:
            continue
        idx = tet.index(0)
        tri = tet[:idx]+ tet[(idx+1):]
        assert 3 == len(tri)
        triangles.append(tri)
    return triangles

"""
May not be properly oriented
"""
def generate_infinite_tetrahedra(mirror):
    # Add all triangles that are not associated with
    # two tets
    infinite_tetrahedra = []
    for (key,value) in mirror.items():
        if 1 == len(value):
            tet = [0] + list(key) # 0 is the "node at infinity"
            infinite_tetrahedra.append(tet)
        else:
            assert 2 == len(value)
    return infinite_tetrahedra

def remove_infinite_tetrahedra(tetrahedra):
    return [tet for tet in tetrahedra if 0 not in tet]

##################################################################################
# Command line stuff

def convert_ctri_to_mesh(ctrifile,meshfile):
    print 'Reading in from CGAL::Triangulation_3 format...'
    (vertices,tetrahedra,neighbors) = read_cgal_tri(ctrifile)
    print '\tNumber of vertices:',len(vertices)-1
    print '\tTotal number of tetrahedra:',len(tetrahedra)
    triangles = generate_boundary_triangles(vertices,tetrahedra)
    print '\tNumber of implied triangles:',len(triangles)
    tetrahedra = remove_infinite_tetrahedra(tetrahedra)
    print '\tNumber of finite tetrahedra:',len(tetrahedra)
    
    # Shouldn't be necessary
    mirror = build_mirror(tetrahedra)
    reorient_tetrahedras(vertices,tetrahedra,mirror)

    ext = meshfile.rsplit('.',1)[1]
    assert('mesh' == ext)
    
    write_inria_mesh(meshfile,vertices,triangles,tetrahedra)

def convert_mesh_to_ctri(meshfile,ctrifile):
    print 'Reading in from INRIA .mesh (Medit) format...'
    (vertices,triangles,tetrahedra) = read_inria_mesh(meshfile)
    print '\tNumber of vertices:',len(vertices)-1
    print '\tNumber of triangles:',len(triangles)
    print '\tNumber of tetrahedra:',len(tetrahedra)
    mirror = build_mirror(tetrahedra)
    infinite_tetrahedra = generate_infinite_tetrahedra(mirror)
    tetrahedra.extend(infinite_tetrahedra)
    print '\tTotal number of tetrahedra:',len(tetrahedra)

    # Rebuild the mirror
    mirror = build_mirror(tetrahedra)
    reorient_tetrahedras(vertices,tetrahedra,mirror)
    neighbors = generate_neighbors(tetrahedra,mirror)

    ext = ctrifile.rsplit('.',1)[1]
    assert('ctri' == ext)
    
    write_cgal_tri(ctrifile,vertices,tetrahedra,neighbors)      

if __name__ == "__main__":

    (_,infile,outfile) = sys.argv

    # Read in files and generate missing information
    in_ext = infile.rsplit('.',1)[1]
    if in_ext == 'ctri':
        convert_ctri_to_mesh(infile,outfile)
    elif in_ext == 'mesh':
        convert_mesh_to_ctri(infile,outfile)
    else:
        print "Unrecognized input extension:",in_ext
        quit()