-
-
Save acrosby/11180502 to your computer and use it in GitHub Desktop.
import numpy as np | |
import struct | |
def read_stl(filename): | |
with open(filename) as f: | |
Header = f.read(80) | |
nn = f.read(4) | |
Numtri = struct.unpack('i', nn)[0] | |
record_dtype = np.dtype([ | |
('Normals', np.float32, (3,)), | |
('Vertex1', np.float32, (3,)), | |
('Vertex2', np.float32, (3,)), | |
('Vertex3', np.float32, (3,)), | |
('atttr', '<i2', (1,)), | |
]) | |
data = np.zeros((Numtri,), dtype=record_dtype) | |
for i in range(0, Numtri, 10): | |
d = np.fromfile(f, dtype=record_dtype, count=10) | |
data[i:i+len(d)] = d | |
#normals = data['Normals'] | |
v1 = data['Vertex1'] | |
v2 = data['Vertex2'] | |
v3 = data['Vertex3'] | |
points = np.hstack(((v1[:, np.newaxis, :]), (v2[:, np.newaxis, :]), (v3[:, np.newaxis, :]))) | |
return points | |
def calc_volume(points): | |
''' | |
Calculate the volume of an stl represented in m x 3 x 3 points array. Expected that input units is mm, so that | |
output is in cubic centimeters (cc). | |
''' | |
v = points | |
volume = np.asarray([np.cross(v[i, 0, :], v[i, 1, :]).dot(v[i, 2, :]) for i in range(points.shape[0])]) | |
return np.abs(volume.sum()/6.)/(10.**3.) | |
def bounding_box(points): | |
''' | |
Calculate the bounding box edge lengths of an stl using the design coordinate system (not an object oriented bounding box), | |
expect that input coordinates are in mm. | |
''' | |
v = points | |
x = v[..., 0].flatten() | |
y = v[..., 1].flatten() | |
z = v[..., 2].flatten() | |
return (x.max()-x.min(), y.max()-y.min(), z.max()-z.min()) | |
def iter_calc_volume(filename): | |
# open as binary | |
with open(filename, 'rb') as f: | |
Header = f.read(80) | |
nn = f.read(4) | |
Numtri = struct.unpack('i', nn)[0] | |
record_dtype = np.dtype([ | |
('Normals', np.float32, (3,)), | |
('Vertex1', np.float32, (3,)), | |
('Vertex2', np.float32, (3,)), | |
('Vertex3', np.float32, (3,)), | |
('atttr', '<i2', (1,)), | |
]) | |
volume = 0. | |
for i in range(0, Numtri, 1): | |
d = np.fromfile(f, dtype=record_dtype, count=1) | |
v1 = d['Vertex1'][0] | |
v2 = d['Vertex2'][0] | |
v3 = d['Vertex3'][0] | |
volume += np.cross(v1, v2).dot(v3) | |
return np.abs(volume/6.)/(10.**3.) | |
def iter_calc_bounding(filename): | |
# open as binary | |
with open(filename, 'rb') as f: | |
Header = f.read(80) | |
nn = f.read(4) | |
Numtri = struct.unpack('i', nn)[0] | |
record_dtype = np.dtype([ | |
('Normals', np.float32, (3,)), | |
('Vertex1', np.float32, (3,)), | |
('Vertex2', np.float32, (3,)), | |
('Vertex3', np.float32, (3,)), | |
('atttr', '<i2', (1,)), | |
]) | |
xmax = -9999 | |
xmin = 9999 | |
ymax = -9999 | |
ymin = 9999 | |
zmax = -9999 | |
zmin = 9999 | |
for i in range(0, Numtri, 1): | |
d = np.fromfile(f, dtype=record_dtype, count=1) | |
v1 = d['Vertex1'] | |
v2 = d['Vertex2'] | |
v3 = d['Vertex3'] | |
v = np.hstack(((v1[:, np.newaxis, :]), (v2[:, np.newaxis, :]), (v3[:, np.newaxis, :]))) | |
x = v[..., 0].flatten() | |
y = v[..., 1].flatten() | |
z = v[..., 2].flatten() | |
tmp_xmin = x.min() | |
tmp_xmax = x.max() | |
tmp_ymin = y.min() | |
tmp_ymax = y.max() | |
tmp_zmin = z.min() | |
tmp_zmax = z.max() | |
xmax = max((tmp_xmax, xmax)) | |
xmin = min((tmp_xmin, xmin)) | |
ymax = max((tmp_ymax, ymax)) | |
ymin = min((tmp_ymin, ymin)) | |
zmax = max((tmp_zmax, zmax)) | |
zmin = min((tmp_zmin, zmin)) | |
X = xmax-xmin | |
Y = ymax-ymin | |
Z = zmax-zmin | |
return (X, Y, Z) | |
def test(filename): | |
points = read_stl(filename) | |
vol1 = calc_volume(points) | |
vol2 = iter_calc_volume(filename) | |
print(vol1, vol2) | |
assert np.isclose(vol1, vol2) # test for equal taking into account | |
# floating point precision | |
bb1 = bounding_box(points) | |
bb2 = iter_calc_bounding(filename) | |
print(bb1, bb2) | |
assert np.isclose(bb1[0], bb2[0]) and np.isclose(bb1[1], bb2[1]) and np.isclose(bb1[2], bb2[2]) | |
return True | |
if __name__ == "__main__": | |
import sys | |
if sys.argv[1] == '-t': | |
test(sys.argv[2]) | |
print("Passed!") | |
else: | |
filename = sys.argv[1] | |
print("Calculating the volume of %s in cc's" % (filename,)) | |
print("The volume is: %f" % (iter_calc_volume(filename),)) | |
print("The bounding box is: %s" % (iter_calc_bounding(filename),)) |
Hi, im trying to write a STL file in Binary format and I cant... can you please help me?
`def surf2stl(filename=None, x=None, y=None, z=None):
p1=[]
p2=[]
p3=[]
if filename=="":
print("error en el nombre de archivo")
fid = open(filename,"w")
titulo_str='Created by Tomas Ferrari'
fid.write("solid %s\n" % titulo_str)
nfacets = 0
a=0.006505
b=-0.0066
c=-0.00391
old_v3=[a ,b, c]
for i in range(0,len(z)-1):
for j in range(0,len(z[0])-1):
p1=[x[i][j],y[i][j],z[i][j]]
p2=[x[i][j+1],y[i][j+1],z[i][j+1]]
p3=[x[i+1][j+1],y[i+1][j+1],z[i+1][j+1]]
valor1=arr_nan(p1,3)
valor2=arr_nan(p2,3)
valor3=arr_nan(p3,3)
valor_tot=valor1+valor2+valor3
if valor_tot!=0:
val=0
else:
val=1
v1=[p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]]
v2=[p3[0]-p1[0],p3[1]-p1[1],p3[2]-p1[2]]
v3 = np.cross(v1, v2)
if v3[0]==0 and v3[1]==0 and v3[2]==0:
v3=old_v3
else:
old_v3=v3
mult = []
total=0
for k in range(len(v3)):
mult.append(v3[k]*v3[k])
total=np.sum(mult)
total=float(total)
tt=math.sqrt(total)
n=[]
for l in range(len(v3)):
n.append(v3[l]/tt)
for m in range(len(n)):
n[m]=float(n[m])
for o in range(len(p1)):
p1[o]=float(p1[o])
for p in range(len(p2)):
p2[p]=float(p2[p])
for q in range(len(p3)):
p3[q]=float(p3[q])
fid.write("facet normal %.7f %.7f %.7f\n" % (n[0], n[1], n[2]))
fid.write("outer loop\n")
fid.write("vertex %.7f %.7f %.7f\n" % (p1[0], p1[1], p1[2]))
fid.write("vertex %.7f %.7f %.7f\n" % (p2[0], p2[1], p2[2]))
fid.write("vertex %.7f %.7f %.7f\n" % (p3[0], p3[1], p3[2]))
fid.write("endloop\n")
fid.write("endfacet\n")
nfacets = nfacets + val
#----------------------------------------------------------------------#
p1=[x[i+1][j+1],y[i+1][j+1],z[i+1][j+1]]
p2=[x[i+1][j],y[i+1][j],z[i+1][j]]
p3=[x[i][j],y[i][j],z[i][j]]
valor1=arr_nan(p1,3)
valor2=arr_nan(p2,3)
valor3=arr_nan(p3,3)
valor_tot=valor1+valor2+valor3
if valor_tot!=0:
val=1
else:
val=0
v1=[p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]]
v2=[p3[0]-p1[0],p3[1]-p1[1],p3[2]-p1[2]]
v3 = np.cross(v1, v2)
if v3[0]==0 and v3[1]==0 and v3[2]==0:
v3=old_v3
else:
old_v3=v3
mult = []
total=0
for k in range(len(v3)):
mult.append(v3[k]*v3[k])
total=np.sum(mult)
total=float(total)
tt=math.sqrt(total)
n=[]
for l in range(len(v3)):
n.append(v3[l]/tt)
for m in range(len(n)):
n[m]=float(n[m])
for o in range(len(p1)):
p1[o]=float(p1[o])
for p in range(len(p2)):
p2[p]=float(p2[p])
for q in range(len(p3)):
p3[q]=float(p3[q])
fid.write("facet normal %.7f %.7f %.7f\n" % (n[0], n[1], n[2]))
fid.write("outer loop\n")
fid.write("vertex %.7f %.7f %.7f\n" % (p1[0], p1[1], p1[2]))
fid.write("vertex %.7f %.7f %.7f\n" % (p2[0], p2[1], p2[2]))
fid.write("vertex %.7f %.7f %.7f\n" % (p3[0], p3[1], p3[2]))
fid.write("endloop\n")
fid.write("endfacet\n")
nfacets = nfacets + val
#----------------------------------------------------------------------#
fid.write("endsolid %s\n" % titulo_str)
print("wrote %d facets" % nfacets)
fid.close()
return`
I wrote this to make ASCII files but need to do it in binary....
Hi, it seems like reading stl file to retrieve the triangular' info including vertices and normal vector was done nicely in your code, but somehow I don't see how i can get the connectivity information of tetrahedron in the stl file. Do you know a technique to retrieve this info? Basically, I need to find out all tetrahedron's connections. For instance, I need to find what triangles made up of a specific tetrahedron in the stl file.
Hi @trtpg, I'm not sure I understand what you are asking. In the case of the STL files that I usually deal with the triangles each represent an element of the outer of surface of the object. Together they are a triangular mesh specifying the surface boundary of a 3d shape. The interior of an object is not represented here, other than being contained by the outer surface.
Hi @acrosby, What I am looking for are connectivities of the mesh, the indices of the vertices describing the tetrahedron. For example if vertices 1, 5, 8, 11 make up a tetrahedron, then I need to know this order and put in a array-like [1,5, 8, 11]. So if there are n tetrahedrons in the mesh then I need a n x 4 list to store this info such as [[1,5, 8, 11], [etc,], .... []]
Hi @acrosby! Your script saved my day, I had two folders full of STLs that I needed to process to find the bounding box.
I've converted this file using the 2to3 tool that comes bundled with Python, then made two very small adjustments to open the files as binary using the "rb" (read binary) mode. That's because Python 3 tries to open the files as text and gives a UnicodeError if it finds something it doesn't like.
@carribeiro I'm glad it was helpful! I haven't used this since the switch to Python3, if you want to post a diff or just or your changed version here I'll update my version of the Gist.
Hi @acrosby here follows my update. I've done the 2to3 upgrade and changed the file mode at the open()
call.
import numpy as np
import struct
def read_stl(filename):
with open(filename) as f:
Header = f.read(80)
nn = f.read(4)
Numtri = struct.unpack('i', nn)[0]
record_dtype = np.dtype([
('Normals', np.float32, (3,)),
('Vertex1', np.float32, (3,)),
('Vertex2', np.float32, (3,)),
('Vertex3', np.float32, (3,)),
('atttr', '<i2', (1,)),
])
data = np.zeros((Numtri,), dtype=record_dtype)
for i in range(0, Numtri, 10):
d = np.fromfile(f, dtype=record_dtype, count=10)
data[i:i+len(d)] = d
#normals = data['Normals']
v1 = data['Vertex1']
v2 = data['Vertex2']
v3 = data['Vertex3']
points = np.hstack(((v1[:, np.newaxis, :]), (v2[:, np.newaxis, :]), (v3[:, np.newaxis, :])))
return points
def calc_volume(points):
'''
Calculate the volume of an stl represented in m x 3 x 3 points array. Expected that input units is mm, so that
output is in cubic centimeters (cc).
'''
v = points
volume = np.asarray([np.cross(v[i, 0, :], v[i, 1, :]).dot(v[i, 2, :]) for i in range(points.shape[0])])
return np.abs(volume.sum()/6.)/(10.**3.)
def bounding_box(points):
'''
Calculate the bounding box edge lengths of an stl using the design coordinate system (not an object oriented bounding box),
expect that input coordinates are in mm.
'''
v = points
x = v[..., 0].flatten()
y = v[..., 1].flatten()
z = v[..., 2].flatten()
return (x.max()-x.min(), y.max()-y.min(), z.max()-z.min())
def iter_calc_volume(filename):
# open as binary
with open(filename, 'rb') as f:
Header = f.read(80)
nn = f.read(4)
Numtri = struct.unpack('i', nn)[0]
record_dtype = np.dtype([
('Normals', np.float32, (3,)),
('Vertex1', np.float32, (3,)),
('Vertex2', np.float32, (3,)),
('Vertex3', np.float32, (3,)),
('atttr', '<i2', (1,)),
])
volume = 0.
for i in range(0, Numtri, 1):
d = np.fromfile(f, dtype=record_dtype, count=1)
v1 = d['Vertex1'][0]
v2 = d['Vertex2'][0]
v3 = d['Vertex3'][0]
volume += np.cross(v1, v2).dot(v3)
return np.abs(volume/6.)/(10.**3.)
def iter_calc_bounding(filename):
# open as binary
with open(filename, 'rb') as f:
Header = f.read(80)
nn = f.read(4)
Numtri = struct.unpack('i', nn)[0]
record_dtype = np.dtype([
('Normals', np.float32, (3,)),
('Vertex1', np.float32, (3,)),
('Vertex2', np.float32, (3,)),
('Vertex3', np.float32, (3,)),
('atttr', '<i2', (1,)),
])
xmax = -9999
xmin = 9999
ymax = -9999
ymin = 9999
zmax = -9999
zmin = 9999
for i in range(0, Numtri, 1):
d = np.fromfile(f, dtype=record_dtype, count=1)
v1 = d['Vertex1']
v2 = d['Vertex2']
v3 = d['Vertex3']
v = np.hstack(((v1[:, np.newaxis, :]), (v2[:, np.newaxis, :]), (v3[:, np.newaxis, :])))
x = v[..., 0].flatten()
y = v[..., 1].flatten()
z = v[..., 2].flatten()
tmp_xmin = x.min()
tmp_xmax = x.max()
tmp_ymin = y.min()
tmp_ymax = y.max()
tmp_zmin = z.min()
tmp_zmax = z.max()
xmax = max((tmp_xmax, xmax))
xmin = min((tmp_xmin, xmin))
ymax = max((tmp_ymax, ymax))
ymin = min((tmp_ymin, ymin))
zmax = max((tmp_zmax, zmax))
zmin = min((tmp_zmin, zmin))
X = xmax-xmin
Y = ymax-ymin
Z = zmax-zmin
return (X, Y, Z)
def test(filename):
points = read_stl(filename)
vol1 = calc_volume(points)
vol2 = iter_calc_volume(filename)
print(vol1, vol2)
assert np.isclose(vol1, vol2) # test for equal taking into account
# floating point precision
bb1 = bounding_box(points)
bb2 = iter_calc_bounding(filename)
print(bb1, bb2)
assert np.isclose(bb1[0], bb2[0]) and np.isclose(bb1[1], bb2[1]) and np.isclose(bb1[2], bb2[2])
return True
if __name__ == "__main__":
import sys
if sys.argv[1] == '-t':
test(sys.argv[2])
print("Passed!")
else:
filename = sys.argv[1]
print("Calculating the volume of %s in cc's" % (filename,))
print("The volume is: %f" % (iter_calc_volume(filename),))
print("The bounding box is: %s" % (iter_calc_bounding(filename),))
Can anyone help me? I get an error:
File "D:_CastingSystem\calc_volume.py", line 7, in read_stl
Header = f.read(80)
^^^^^^^^^^
File "C:\Users\New\AppData\Local\Programs\Python\Python311\Lib\encodings\cp1251.py", line 23, in decode
return codecs.charmap_decode(input,self.errors,decoding_table)[0]
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
UnicodeDecodeError: 'charmap' codec can't decode byte 0x98 in position 152: character maps to
Thanks, I forgot this existed!