1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
| from __future__ import with_statement
import numpy as np
def loadSTL(filename):
"""Parse and load and ASCII STL file.
Returns a list of triangles and a list of normals."""
# could use some error checking
triangles = []
normals = []
t = []
with open(filename) as f:
for line in f:
line = line.lstrip()
if line.startswith('facet normal'):
normals.append(map(float, line.split()[2:]))
elif line.startswith('vertex'):
t.append(map(float, line.split()[1:]))
if len(t) == 3:
triangles.append(t)
t = []
return triangles, normals
def stl2zmap(filename,(xstep,ystep)=(1.0,1.0)):
"""Returns a height map relative to the z=0 plane"""
t,n = loadSTL(model)
t = np.array(t)
n = np.array(n)
# remove triangles facing away
b = n[:,2] > 0 # normal.z > 0
n = n.compress(b,0)
t = t.compress(b,0)
# calculate bounds
xmax, ymax, zmax = np.max(np.max(t,0),0)
xmin, ymin, zmin = np.min(np.min(t,0),0)
# calculate 2D vectors for each triangle
# used for triangle inclusion test
A,B,C = t.transpose(1,0,2)[:,:,:2]
AB = (B-A)
BC = (C-B)
CA = (A-C)
# divide normals by their z component
# used to compute z values
nz = n[:,2]
nz = np.array([nz,nz,nz]).T
n /= nz
# main loop
data = []
for y in np.arange(ymin, ymax+ystep, ystep):
line = []
data.append(line)
for x in np.arange(xmin, xmax+xstep, xstep):
# triangle inclusion test: check that the point
# is on the same side of each vector of the triangle
s1 = np.cross([x,y] - A, AB) >= 0
s2 = np.cross([x,y] - B, BC) >= 0
s3 = np.cross([x,y] - C, CA) >= 0
b = (s1 == s2) & (s2 == s3)
# keep only these triangles
if b.any():
# we just need one point from each triangle
pt = t.compress(b,0).transpose(1,0,2)[0]
# and normals
nn = n.compress(b,0)
# calculate z for each remaining triangle
z = np.sum(nn * (pt - [x,y,0]),1) # z = (nx/nz)(x0-x)+(ny/nz)(y0-y)+z0
# get maximum value
line.append(np.max(z))
else:
line.append(zmin)
return data
def zmap2img(data):
"""Convert a 2D floating point list to a greyscale image."""
import Image
data = np.array(data)
# initialize image
ysize, xsize = data.shape
im = Image.new('L',(xsize, ysize))
pix = im.load()
# get bounds
zmin, zmax = np.min(data), np.max(data)
scale = 255 / (zmax-zmin)
# scale and write to image
for y,line in enumerate(data):
for x,z in enumerate(line):
pix[x,ysize-y-1] = (z - zmin) * scale + 0.5
return im
if __name__ == '__main__':
model = "C:\\Prog\\python\\3D\\bottle.stl"
data = stl2zmap(model,(0.2,0.2))
# make a PIL Image
img = zmap2img(data)
img.show() |
