diff --git a/Scratch.ipynb b/Scratch.ipynb index 4d91923..719be43 100644 --- a/Scratch.ipynb +++ b/Scratch.ipynb @@ -6,22 +6,36 @@ "source": [ "To do:\n", "\n", - "- Parametrize gen_twisted_boundary over boundaries\n", + "- Parametrize gen_twisted_boundary over boundaries and\n", + "do my nested spiral\n", + "- Rewrite ram_horn in terms of newer abstractions\n", "- Encode the notions of \"generator which transforms an\n", "existing list of boundaries\", \"generator which transforms\n", "another generator\"\n", "- Work directly with lists of boundaries. The only thing\n", "I ever do with them is apply transforms to all of them, or\n", "join adjacent ones with corresponding elements.\n", - "- Refactor/rewrite 'concat' with meshes because I think it\n", - "is scaling horribly. Perhaps have pre-allocated arrays that\n", - "scale by a factor of 1.5 every time they need to be grown, or\n", - "something.\n", "- Why do I get the weird zig-zag pattern on the triangles,\n", "despite larger numbers of them? Is it something in how I\n", - "twist the frames? (Note that I always connect one diagonal\n", - "on the boundaries - but I can choose the other. This might\n", - "not matter, but I should understand why.)" + "twist the frames?\n", + " - How can I compute the *torsion* on a quad? I think it\n", + " comes down to this: torsion applied across the quad I'm\n", + " triangulating leading to neither diagonal being a\n", + " particularly good choice. Subdividing the boundary seems\n", + " to help, but other triangulation methods (e.g. turning a\n", + " quad to 4 triangles by adding the centroid) could be good\n", + " too.\n", + " - Facets/edges are just oriented the wrong way...\n", + "- I need an actual example of branching/forking. If I simply\n", + "split a boundary into sub-boundaries per the rules I already\n", + "have in my notes, then this still lets me split any way I want\n", + "to without having to worry about joining N boundaries instead\n", + "of 2, doesn't it?\n", + "\n", + "Other notes:\n", + "- Picking at random the diagonal on the quad to triangulate with\n", + " does seem to turn 'error' just to noise, and in its own way this\n", + " is preferable. " ] }, { @@ -41,7 +55,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -49,24 +63,24 @@ "#mesh = examples.twist()\n", "#mesh = examples.twist_nonlinear()\n", "#mesh = examples.twist_from_gen()\n", - "#mesh = examples.twisty_torus(4000)" + "mesh = examples.twisty_torus(1000)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# for me to be able to nest this, gen_twisted_boundary has to\n", "# be parametrized over some other boundary (the way )\n", - "gen = examples.gen_inc_y(examples.gen_twisted_boundary())\n", - "mesh = gen2mesh(gen, 100, True)" + "#gen = examples.gen_inc_y(examples.gen_twisted_boundary())\n", + "#mesh = gen2mesh(gen, 100, True)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -75,8 +89,9 @@ "center = meshutil.Transform().translate(-0.5, -0.5, -0.5)\n", "base = meshutil.cube(open_xz=False).transform(center)\n", "base = base.transform(wtf)\n", + "mesh_out = mesh\n", "# to enable:\n", - "mesh_out = mesh.concat(base)\n", + "mesh_out = mesh_out.concat(base)\n", "\n", "mesh_fname = \"twist_wtf.stl\"\n", "mesh_out.to_stl_mesh().save(mesh_fname)" @@ -84,9 +99,229 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "face_normals didn't match triangles, ignoring!\n", + "/home/hodapp/.local/lib/python3.6/site-packages/IPython/core/display.py:694: UserWarning: Consider using IPython.display.IFrame instead\n", + " warnings.warn(\"Consider using IPython.display.IFrame instead\")\n" + ] + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "# TODO: Just use the mesh data directly (no sense in saving & re-loading)\n", "m = trimesh.load_mesh(mesh_fname)\n", @@ -94,13 +329,6 @@ "scene = trimesh.Scene([m])\n", "scene.show()" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/examples.py b/examples.py index 922a060..3097516 100755 --- a/examples.py +++ b/examples.py @@ -19,8 +19,9 @@ def ram_horn(): ], dtype=numpy.float64) - [0.5, 0.5, 0] xf0_to_1 = meshutil.Transform().translate(0,0,1) b1 = xf0_to_1.apply_to(b0) - mesh = meshutil.join_boundary_simple(b0, b1) - mesh = mesh.concat(meshutil.close_boundary_simple(b0)) + meshes = [] + meshes.append(meshutil.join_boundary_simple(b0, b1)) + meshes.append(meshutil.close_boundary_simple(b0)) for i in range(4): # Opening boundary: b = b1 @@ -37,10 +38,10 @@ def ram_horn(): .translate(0,0,0.8) b_sub1 = incr.compose(xf).apply_to(b) m = meshutil.join_boundary_simple(b_sub0, b_sub1) - mesh = mesh.concat(m) + meshes.append(m) xf = incr.compose(xf) # Close final boundary: - mesh = mesh.concat(meshutil.close_boundary_simple(b_sub1[::-1,:])) + meshes.append(meshutil.close_boundary_simple(b_sub1[::-1,:])) # ::-1 is to reverse the boundary's order to fix winding order. # Not sure of the "right" way to fix winding order here. # The boundary vertices go in an identical order... it's just @@ -51,6 +52,7 @@ def ram_horn(): # I don't need to subdivide *geometry*. # I need to subdivide *space* and then put geometry in it. + mesh = meshutil.FaceVertexMesh.concat_many(meshes) return mesh # Interlocking twists. @@ -63,17 +65,13 @@ def twist(ang=0.1, dz=0.2, dx0=2, count=4, scale=0.98): [1, 1, 0], [0, 1, 0], ], dtype=numpy.float64) - [0.5, 0.5, 0] - mesh = None + meshes = [] for i in range(count): xf = meshutil.Transform() \ .translate(dx0, 0, 0) \ .rotate([0,0,1], numpy.pi * 2 * i / count) b0 = xf.apply_to(b) - m = meshutil.close_boundary_simple(b0) - if mesh is None: - mesh = m - else: - mesh = mesh.concat(m) + meshes.append(meshutil.close_boundary_simple(b0)) for layer in range(256): b_sub0 = xf.apply_to(b) incr = meshutil.Transform() \ @@ -82,10 +80,11 @@ def twist(ang=0.1, dz=0.2, dx0=2, count=4, scale=0.98): .scale(scale) b_sub1 = xf.compose(incr).apply_to(b) m = meshutil.join_boundary_simple(b_sub0, b_sub1) - mesh = mesh.concat(m) + meshes.append(m) xf = xf.compose(incr) # Close final boundary: - mesh = mesh.concat(meshutil.close_boundary_simple(b_sub1[::-1,:])) + meshes.append(meshutil.close_boundary_simple(b_sub1[::-1,:])) + mesh = meshutil.FaceVertexMesh.concat_many(meshes) return mesh def twist_nonlinear(dx0 = 2, dz=0.2, count=3, scale=0.99, layers=100): @@ -100,17 +99,13 @@ def twist_nonlinear(dx0 = 2, dz=0.2, count=3, scale=0.99, layers=100): [1, 1, 0], [0, 1, 0], ], dtype=numpy.float64) - [0.5, 0.5, 0] - mesh = None + meshes = [] for i in range(count): xf = meshutil.Transform() \ .translate(dx0, 0, 0) \ .rotate([0,0,1], numpy.pi * 2 * i / count) b0 = xf.apply_to(b) - m = meshutil.close_boundary_simple(b0) - if mesh is None: - mesh = m - else: - mesh = mesh.concat(m) + meshes.append(meshutil.close_boundary_simple(b0)) for layer in range(layers): b_sub0 = xf.apply_to(b) ang = ang_fn(layer) @@ -120,10 +115,11 @@ def twist_nonlinear(dx0 = 2, dz=0.2, count=3, scale=0.99, layers=100): .scale(scale) b_sub1 = xf.compose(incr).apply_to(b) m = meshutil.join_boundary_simple(b_sub0, b_sub1) - mesh = mesh.concat(m) + meshes.append(m) xf = xf.compose(incr) # Close final boundary: - mesh = mesh.concat(meshutil.close_boundary_simple(b_sub1[::-1,:])) + meshes.append(meshutil.close_boundary_simple(b_sub1[::-1,:])) + mesh = meshutil.FaceVertexMesh.concat_many(meshes) return mesh # Generate a frame with 'count' boundaries in the XZ plane. @@ -138,6 +134,7 @@ def gen_twisted_boundary(count=4, dx0=2, dz=0.2, ang=0.1): ], dtype=numpy.float64) - [0.5, 0, 0.5] b = meshutil.subdivide_boundary(b) b = meshutil.subdivide_boundary(b) + b = meshutil.subdivide_boundary(b) # Generate 'seed' transformations: xfs = [meshutil.Transform().translate(dx0, 0, 0).rotate([0,1,0], numpy.pi * 2 * i / count) for i in range(count)] @@ -183,29 +180,30 @@ def gen_torus_xy(gen, rad=2, frames=100): # String together boundaries from a generator. # If count is nonzero, run only this many iterations. -def gen2mesh(gen, count=0, flip_order=False, loop=False): +def gen2mesh(gen, count=0, flip_order=False, loop=False, join_fn=meshutil.join_boundary_optim): # Get first list of boundaries: bs_first = next(gen) bs_last = bs_first # TODO: Begin and end with close_boundary_simple - mesh = meshutil.FaceVertexMesh.Empty() + meshes = [] for i,bs_cur in enumerate(gen): if count > 0 and i >= count: break for j,b in enumerate(bs_cur): if flip_order: - m = meshutil.join_boundary_simple(b, bs_last[j]) + m = join_fn(b, bs_last[j]) else: - m = meshutil.join_boundary_simple(bs_last[j], b) - mesh = mesh.concat(m) + m = join_fn(bs_last[j], b) + meshes.append(m) bs_last = bs_cur if loop: for b0,b1 in zip(bs_last, bs_first): if flip_order: - m = meshutil.join_boundary_simple(b1, b0) + m = join_fn(b1, b0) else: - m = meshutil.join_boundary_simple(b0, b1) - mesh = mesh.concat(m) + m = join_fn(b0, b1) + meshes.append(m) + mesh = meshutil.FaceVertexMesh.concat_many(meshes) return mesh def twist_from_gen(): @@ -216,12 +214,21 @@ def twist_from_gen(): # frames = How many step to build this from: # turn = How many full turns to make in inner twist # count = How many inner twists to have -def twisty_torus(frames = 200, turns = 4, count = 4, rad = 4): +def twisty_torus(frames = 5000, turns = 4, count = 4, rad = 4): # In order to make this line up properly: angle = numpy.pi * 2 * turns / frames gen = gen_torus_xy(gen_twisted_boundary(count=count, ang=angle), rad=rad, frames=frames) return gen2mesh(gen, 0, flip_order=True, loop=True) +# frames = How many step to build this from: +# turn = How many full turns to make in inner twist +# count = How many inner twists to have +def twisty_torus_opt(frames = 200, turns = 4, count = 4, rad = 4): + # In order to make this line up properly: + angle = numpy.pi * 2 * turns / frames + gen = gen_torus_xy(gen_twisted_boundary(count=count, ang=angle), rad=rad, frames=frames) + return gen2mesh(gen, 0, flip_order=True, loop=True, join_fn=meshutil.join_boundary_optim) + def main(): fns = { ram_horn: "ramhorn.stl", diff --git a/meshutil.py b/meshutil.py index 5dd7e90..0b4e902 100644 --- a/meshutil.py +++ b/meshutil.py @@ -38,9 +38,28 @@ class FaceVertexMesh(object): v[i] = [self.v[iv0], self.v[iv1], self.v[iv2]] return stl.mesh.Mesh(data) @classmethod - def Empty(self): + def Empty(cls): return FaceVertexMesh(numpy.zeros((0,3)), numpy.zeros((0,3), dtype=int)) - + @classmethod + def concat_many(cls, meshes): + nv = 0 + nf = 0 + for m in meshes: + nv += m.v.shape[0] + nf += m.f.shape[0] + v = numpy.zeros((nv,3), dtype=numpy.float64) + f = numpy.zeros((nf,3), dtype=int) + vi = 0 + fi = 0 + for m in meshes: + vj = vi + m.v.shape[0] + fj = fi + m.f.shape[0] + v[vi:vj,:] = m.v + f[fi:fj,:] = m.f + vi + vi = vj + fi = fj + return FaceVertexMesh(v, f) + class Transform(object): def __init__(self, mtx=None): if mtx is None: @@ -207,6 +226,16 @@ def join_boundary_simple(bound1, bound2): fs[2*i + 1] = [v1, n + v0, v0] return FaceVertexMesh(vs, fs) +def join_boundary_optim(bound1, bound2): + # bound1 and bound2 must stay in order, but we can rotate + # the starting point to whatever we want. Use distance as + # a metric: + errs = [numpy.linalg.norm(bound1 - numpy.roll(bound2, i, axis=0)) + for i,_ in enumerate(bound1)] + # What shift gives the lowest distance? + i = numpy.argmin(errs) + return join_boundary_simple(bound1, numpy.roll(bound2, i, axis=0)) + def close_boundary_simple(bound): # This will fail for any non-convex boundary! centroid = numpy.mean(bound, axis=0)