git.haldean.org symrep / 9c6b7e2
more geometry tools Will Haldean Brown 6 years ago
3 changed file(s) with 64 addition(s) and 24 deletion(s). Raw diff Collapse all Expand all
00 import numpy as np
11 from symrep import *
2 from symrep.solids import *
23
3 n = solids.union(
4 solids.sphere(const(1)),
5 solids.intersection(
6 solids.sphere(const(3)),
7 solids.translate(
8 solids.sphere(const(2)),
9 const(np.array((3., 0., 0.))),
4 n = union(
5 sphere(const(1)),
6 intersection(
7 sphere(const(3)),
8 translate(
9 sphere(const(2)),
10 vec(3., 0., 0.),
1011 ),
12 ),
13 translate(
14 difference(
15 sphere(const(2)),
16 sphere(const(1.5)),
17 translate(
18 sphere(const(2)),
19 vec(1., 0., 0.),
20 ),
21 ),
22 vec(0., 3., 0.),
1123 ),
1224 )
1325
26 with open("sphere.dot", "w") as f:
27 to_dot(n, f)
28
29 points = sample_solid(
30 n, (-2, -2, -2), (5, 5, 2), 20000, is_inside)
1431 with open("test.xyz", "w") as f:
15 solids.write_pointcloud(
16 n, (-1, -2, -2), (5, 2, 2), 100000, f)
32 solids.write_point_cloud(points, f)
00 import functools
1 import numpy as np
12 import operator
3 import random
24
35
46 def const(val):
11 import random
22 import symrep.base
33
4
5 def vec(*v):
6 if len(v) == 3:
7 fill = np.zeros(4)
8 fill[:3] = v
9 v = fill
10 return symrep.base.const(v)
11
12 def point(*p):
13 if len(p) == 3:
14 fill = np.ones(4)
15 fill[:3] = p
16 p = fill
17 return symrep.base.const(p)
418
519 def transform(n, T, name=None, nodetype="transform"):
620 return symrep.base.Node(
721 name, nodetype, lambda t: n(T(t).dot(t)), {n: "value", T: "transform"})
822
923 def translate(n, trans, name=None):
10 def getT(t):
11 T = np.eye(4)
12 T[:3, 3] = -trans(t)
13 return T
14 return transform(n, getT, name=name, nodetype="translate")
15
16 def sphere(radius, name=None):
17 return symrep.base.Node(
18 name, "sphere",
19 lambda t: np.linalg.norm(t[:3]) - radius(t), {radius: "radius"})
24 return symrep.base.shift(n, invert(trans), name=name)
2025
2126 def union(*nodes):
2227 return symrep.base.Node(
2833 None, "intersection", lambda t: max(map(lambda n: n(t), nodes)),
2934 {n: "value" for n in nodes})
3035
31 def to_pointcloud(root, hi, lo, num_points):
36 def invert(n, name=None):
37 return symrep.base.Node(name, "invert", lambda t: -n(t), {n: "value"})
38
39 def difference(base, *subs):
40 return intersection(base, invert(union(*subs)))
41
42 def sphere(radius, name=None):
43 return symrep.base.Node(
44 name, "sphere",
45 lambda t: np.linalg.norm(t[:3]) - radius(t), {radius: "radius"})
46
47 def sample_solid(root, lo, hi, num_points, eval_pt):
3248 found = 0
3349 tried = 0
3450 t = np.ones(4)
3652 t[0] = random.uniform(lo[0], hi[0])
3753 t[1] = random.uniform(lo[1], hi[1])
3854 t[2] = random.uniform(lo[2], hi[2])
39 if root(t) < 0:
55 if eval_pt(root, t):
4056 yield t
4157 found += 1
4258 tried += 1
4359 print("tried {} points to find {}".format(tried, num_points))
4460
45 def write_pointcloud(root, hi, lo, num_points, f):
46 num_written = 0
47 for point in to_pointcloud(root, hi, lo, num_points):
61 def is_inside(root, t):
62 return root(t) <= 0
63
64 def is_on_surface(root, t, accuracy):
65 val = root(t)
66 return val <= 0 and abs(val) <= accuracy
67
68 def write_point_cloud(points, f):
69 for point in points:
4870 f.write("{} {} {}\n".format(*(point[:3])))