# How to extract surface triangles from a tetrahedral mesh?

I want to render a tetrahedral mesh using some 3D software. However, I cannot directly load the tetrahedral mesh in the software of my choice (e.g. Blender) as the file format that I have for tetrahedral meshes is not supported. So I should somehow extract the faces with corresponding vertex indices myself.

For a cube, my tetrahedral file contains vertex IDs for each tetrahedron which includes 4 faces is as follow:

```v 0.41 0.41 0.41
v 0.41 0.41 -0.41
v 0.41 -0.41 0.41
v 0.41 -0.41 -0.41
v -0.41 0.41 0.41
v -0.41 0.41 -0.41
v -0.41 -0.41 0.41
v -0.41 -0.41 -0.41
t 0 1 2 4
t 5 1 4 7
t 1 2 4 7
t 3 1 7 2
t 6 4 2 7
```

However, I’m not sure how I can and extract the surface mesh given this data. Does someone know how I can do this or what the algorithm is?

here is a simplistic brute force method. For each tetrahedron, for example look at the third one, `t: 1 2 4 7`, by removing each vertex, generate all four combination of three vertices out of the four tetrahedral vertices, i.e.

```face[t][0]: 1 2 4,  face[t][1]: 1 2 7,  face[t][2]: 1 4 7, face[t][3]: 2 4 7
```

and sort each triangle’s integer labels in ascending order (for uniqueness) This way, you can generate the list (or some kind of array) of all faces of all tetrahedral from the tetrahedral mesh.

Now run a loop over the list of all triangle faces that you have just generated, looking for duplicates. Whenever a triangle is contained twice in the list of all triangle faces, you remove it, because it is an interior triangle, i.e. two adjacent tetrahedral share this triangular face, so it is interior face and not a boundary one.

Whatever is left after this procedure, are only the boundary (i.e. the surface) triangle faces of the tetrahedral mesh.

Here is an example of this algorithm written in python

```import numpy as np

def list_faces(t):
t.sort(axis=1)
n_t, m_t= t.shape
f = np.empty((4*n_t, 3) , dtype=int)
i = 0
for j in range(4):
f[i:i+n_t,0:j] = t[:,0:j]
f[i:i+n_t,j:3] = t[:,j+1:4]
i=i+n_t
return f

def extract_unique_triangles(t):
_, indxs, count  = np.unique(t, axis=0, return_index=True, return_counts=True)
return t[indxs[count==1]]

def extract_surface(t):
f=list_faces(t)
f=extract_unique_triangles(f)
return f

V = np.array([
[ 0.41,  0.41,  0.41],
[ 0.41,  0.41, -0.41],
[ 0.41, -0.41,  0.41],
[ 0.41, -0.41, -0.41],
[-0.41,  0.41,  0.41],
[-0.41,  0.41, -0.41],
[-0.41, -0.41,  0.41],
[-0.41, -0.41, -0.41]])

T = np.array([
[0, 1, 2, 4],
[5, 1, 4, 7],
[1, 2, 4, 7],
[3, 1, 7, 2],
[6, 4, 2, 7]])

F_all = list_faces(T)
print(F_all)
print(F_all.shape)

F_surf = extract_surface(T)
print(F_surf)
print(F_surf.shape)
```