Selection of Face of a STL by Face Normal value Threshold

I want to write a script in Python which can generate facegroups in a STL as per the Face Normal value condition. For example, Provided is the snap of Stl, Different colour signifies the face group containing the triangular faces satisfying my given face normal threshold. Is there any simple way to do this in python? Face Group STL

I'm sure there's a python library to load stl files, but I've always just written my own, since the file format is pretty simple (see the Wikipedia article for file format description).

Here is my code to read the stl file:

import numpy as np
import struct

def Unique(inputList):
      """ 
      Given an M x N list, this function gets the unique rows by treating all
      M Ntuples as single objects. This function also returns the indexing
      to convert the unique returned list back to the original non-unique list.
      """

      hashTable=dict()

      indexList=[]
      uniqueList=[]

      indx=0
      for ntuple in inputList:
            if not ntuple in hashTable:
                hashTable[ntuple]=indx
                indexList.append(indx)
                uniqueList.append(ntuple)
                indx+=1
            else:
                indexList.append(hashTable.get(ntuple))      

      return uniqueList, indexList


def IsBinarySTL(filename):
    try:
        with open(filename,'r') as f:
              test=f.readline()
    except UnicodeDecodeError:
        return True

    if len(test) < 5:
        return True
    elif test[0:5].lower() == 'solid':
        return False  # ASCII STL
    else:
        return True

def ReadSTL(filename):
    """ Returns numpy arrays for vertices and facet indexing """
    def GetListFromASCII(filename):
        """ Returns vertex listing from ASCII STL file """
        outputList=[]

        with open(filename,'r') as f:
            lines=[line.split() for line in f.readlines()]
        for line in lines:
            if line[0] == 'vertex':
                    outputList.append(tuple([float(x) for x in line[1:]]))
        return outputList

    def GetListFromBinary(filename):
        """ Returns vertex listing from binary STL file """
        outputList=[]
        with open(filename,'rb') as f:
            f.seek(80) # skip header
            nFacets=struct.unpack('I',f.read(4))[0] # number of facets in piece

            for i in range(nFacets):
                  f.seek(12,1) # skip normal
                  outputList.append(struct.unpack('fff',f.read(12))) # append each vertex triple to list (each facet has 3 vertices)
                  outputList.append(struct.unpack('fff',f.read(12))) 
                  outputList.append(struct.unpack('fff',f.read(12)))
                  f.seek(2,1) # skip attribute
        return outputList

    if IsBinarySTL(filename):
        vertexList = GetListFromBinary(filename)
    else:
        vertexList = GetListFromASCII(filename)

    coords, tempindxs = Unique(vertexList)

    indxs = list()
    templist = list()
    for i in range(len(tempindxs)):
        if (i > 0 ) and not (i % 3):
            indxs.append(templist)
            templist = list()
        templist.append(tempindxs[i])
    indxs.append(templist)

    return np.array(coords), np.array(indxs)

And here is code to compute the facet normals (assuming right-hand-rule)

def GetNormals(vertices, facets):
    """ Returns normals for each facet of mesh """
    u = vertices[facets[:,1],:] - vertices[facets[:,0],:]
    v = vertices[facets[:,2],:] - vertices[facets[:,0],:]
    normals = np.cross(u,v)
    norms = np.sqrt(np.sum(normals*normals, axis=1))
    return normals/norms[:, np.newaxis]

Finally, code to write out the stl file (assuming a list of attributes for each facet):

def WriteSTL(filename, vertices, facets, attributes, header):
    """
    Writes vertices and facets to an stl file. Notes:
    1.) header can not be longer than 80 characters
    2.) length of attributes must be equal to length of facets
    3.) attributes must be integers
    """
    nspaces = 80 - len(header)
    header += nspaces*'\0'

    nFacets = np.shape(facets)[0]
    stl = vertices[facets,:].tolist()

    with open(filename,'wb') as f: # binary
        f.write(struct.pack('80s', header.encode('utf-8'))) # header
        f.write(struct.pack('I',nFacets)) # number of facets
        for i in range(nFacets):
            f.write(struct.pack('fff',0,0,0)) # normals set to 0
            for j in range(3):
                f.write(struct.pack('fff',stl[i][j][0], stl[i][j][1], stl[i][j][2])) # 3 vertices per facet 
            f.write(struct.pack("H", attributes[i])) # 2-byte attribute

Putting this all together, you can do something like the following:

if __name__ == "__main__":
    filename = "bunny.stl"

    vertices, facets = ReadSTL(filename)  # parse stl file
    normals = GetNormals(vertices, facets)  # compute normals

    # Get some value related to normals
    attributes = []
    for i in range(np.shape(normals)[0]):
        attributes.append(int(255*np.sum(normals[i])**2))

    # Write new stl file
    WriteSTL("output.stl", vertices, facets, attributes, "stlheader")

this code snippet reads an stl file, computes the normals, and then assigns an attribute value based on the squared-sum of each normal (note that the attribute must be an integer).

The input and output of this script look like the following: