Compare commits
No commits in common. "521ad5d8a5c310b4a6deb50a32472042f45e7d01" and "2dc86d23bf4e2f4eb286ddbe23597f5daa7bba6a" have entirely different histories.
521ad5d8a5
...
2dc86d23bf
3
.gitignore
vendored
3
.gitignore
vendored
@ -1,5 +1,4 @@
|
||||
*~
|
||||
#*#
|
||||
__pycache__/
|
||||
.ipynb_checkpoints/
|
||||
.ccls-cache/
|
||||
.ipynb_checkpoints/*
|
||||
|
||||
21
README.md
21
README.md
@ -3,14 +3,6 @@
|
||||
This is repo has a few projects that are related in terms of
|
||||
high-level goal, but almost completely unrelated in their descent.
|
||||
|
||||
- `python_isosurfaces_2018_2019` is some Python & Maxima code from
|
||||
2018-2019 from me trying to turn the usual spiral isosurface into a
|
||||
parametric formula of sorts in order to triangulate it more
|
||||
effectively.
|
||||
- `parallel_transport` is some Python code from 2019 September which
|
||||
implemented parallel frame transport, i.e.
|
||||
[Parallel Transport Approach to Curve Framing](https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.8103). It is mostly scratch
|
||||
code / proof-of-concept.
|
||||
- `python_extrude_meshgen` is some Python code from around 2019
|
||||
September which did a sort of extrusion-based code generation.
|
||||
While this had some good results and some good ideas, the basic
|
||||
@ -32,21 +24,10 @@ high-level goal, but almost completely unrelated in their descent.
|
||||
and found a section I'd ignored on the difficulties of producing
|
||||
good meshes from isosurfaces for the sake of rendering. I kept
|
||||
the code around because I figured it would be useful to refer to
|
||||
later, particularly for the integration with Blender - but
|
||||
otherwise shelved this effort.
|
||||
later, particularly for the integration with Blender.
|
||||
- `blender_scraps` contains some scraps of Python code meant to be
|
||||
used inside of Blender's Python scripting - and it contains some
|
||||
conversions from another project, Prosha, for procedural mesh
|
||||
generation in Rust (itself based on learnings from
|
||||
`python_extrude_meshgen`). These examples were proof-of-concept of
|
||||
generating meshes as control cages rather than as "final" meshes.
|
||||
|
||||
It would probably make sense to rename this repo to something with
|
||||
`procedural` in the name rather than `automata` since at some point it
|
||||
ceased to have much to do with automata.
|
||||
|
||||
## Projects not covered here
|
||||
|
||||
- curl-noise work (both in Clojure and in Python/vispy)
|
||||
- parallel transport
|
||||
- prosha, of course
|
||||
|
||||
22
implicit.org
22
implicit.org
@ -1,22 +0,0 @@
|
||||
* Slush-pile from implicit rant
|
||||
(TODO: This was shelved from another post. Do something with it.)
|
||||
|
||||
None of my attempts have ever quite worked as soon as I start
|
||||
dealing with implicit surfaces that have much sharpness to
|
||||
them... which isn't saying much, since I am not exactly world-class
|
||||
at dealing with meshes and mesh topology directly, but it's also
|
||||
legitimately a very hard problem. Even when I tried in CGAL, which
|
||||
can handle implicit surfaces and does a far better job at generating
|
||||
"good" meshes than anything I could hope to write myself, the
|
||||
problem I ran into is that it was generating good meshes for
|
||||
industrial CAD, not good meshes for /rendering/. That is, their
|
||||
inaccuracy to the true surface was bounded, and I could always
|
||||
reduce it by just using denser meshes. For CAD applications, this
|
||||
seems fine. For rendering, though, I kept running into some major
|
||||
inefficiencies in terms of number of triangles required (and as well
|
||||
the amount of time required). I don't know much about the algorithms
|
||||
they use, but I am fairly sure that for detail level N (inversely
|
||||
proportional to minimum feature size) both triangle count, time, and
|
||||
likely memory usage grow with ~O(N^3)~.
|
||||
|
||||
(TODO: Give an example of the above perhaps?)
|
||||
File diff suppressed because one or more lines are too long
@ -1,84 +0,0 @@
|
||||
import numpy
|
||||
import sympy
|
||||
from sympy.vector import CoordSys3D
|
||||
|
||||
def mtx_rotate_by_vector(b, theta):
|
||||
"""Returns 3x3 matrix for rotating around some 3D vector."""
|
||||
# Source:
|
||||
# Appendix A of "Parallel Transport to Curve Framing"
|
||||
c = numpy.cos(theta)
|
||||
s = numpy.sin(theta)
|
||||
rot = numpy.array([
|
||||
[c+(b[0]**2)*(1-c), b[0]*b[1]*(1-c)-s*b[2], b[2]*b[0]*(1-c)+s*b[1]],
|
||||
[b[0]*b[1]*(1-c)+s*b[2], c+(b[1]**2)*(1-c), b[2]*b[1]*(1-c)-s*b[0]],
|
||||
[b[0]*b[2]*(1-c)-s*b[1], b[1]*b[2]*(1-c)+s*b[0], c+(b[2]**2)*(1-c)],
|
||||
])
|
||||
return rot
|
||||
|
||||
def gen_faces(v1a, v1b, v2a, v2b):
|
||||
"""Returns faces (as arrays of vertices) connecting two pairs of
|
||||
vertices."""
|
||||
# Keep winding order consistent!
|
||||
f1 = numpy.array([v1b, v1a, v2a])
|
||||
f2 = numpy.array([v2b, v1b, v2a])
|
||||
return f1, f2
|
||||
|
||||
def approx_tangents(points):
|
||||
"""Returns an array of approximate tangent vectors. Assumes a
|
||||
closed path, and approximates point I using neighbors I-1 and I+1 -
|
||||
that is, treating the three points as a circle.
|
||||
|
||||
Input:
|
||||
points -- Array of shape (N,3). points[0,:] is assumed to wrap around
|
||||
to points[-1,:].
|
||||
|
||||
Output:
|
||||
tangents -- Array of same shape as 'points'. Each row is normalized.
|
||||
"""
|
||||
d = numpy.roll(points, -1, axis=0) - numpy.roll(points, +1, axis=0)
|
||||
d = d/numpy.linalg.norm(d, axis=1)[:,numpy.newaxis]
|
||||
return d
|
||||
|
||||
def approx_arc_length(points):
|
||||
p2 = numpy.roll(points, -1, axis=0)
|
||||
return numpy.sum(numpy.linalg.norm(points - p2, axis=1))
|
||||
|
||||
def torsion(v, arg):
|
||||
"""Returns an analytical SymPy expression for torsion of a 3D curve.
|
||||
|
||||
Inputs:
|
||||
v -- SymPy expression returning a 3D vector
|
||||
arg -- SymPy symbol for v's variable
|
||||
"""
|
||||
# https://en.wikipedia.org/wiki/Torsion_of_a_curve#Alternative_description
|
||||
dv1 = v.diff(arg)
|
||||
dv2 = dv1.diff(arg)
|
||||
dv3 = dv2.diff(arg)
|
||||
v1_x_v2 = dv1.cross(dv2)
|
||||
# This calls for the square of the norm in denominator - but that
|
||||
# is just dot product with itself:
|
||||
return v1_x_v2.dot(dv3) / (v1_x_v2.dot(v1_x_v2))
|
||||
|
||||
def torsion_integral(curve_expr, var, a, b):
|
||||
# The line integral from section 3.1 of "Parallel Transport to Curve
|
||||
# Framing". This should work in theory, but with the functions I've
|
||||
# actually tried, evalf() is ridiculously slow.
|
||||
c = torsion(curve_expr, var)
|
||||
return sympy.Integral(c * (sympy.diff(curve_expr, var).magnitude()), (var, a, b))
|
||||
|
||||
def torsion_integral_approx(curve_expr, var, a, b, step):
|
||||
# A numerical approximation of the line integral from section 3.1 of
|
||||
# "Parallel Transport to Curve Framing"
|
||||
N = CoordSys3D('N')
|
||||
# Get a (callable) derivative function of the curve:
|
||||
curve_diff = curve_expr.diff(var)
|
||||
diff_fn = sympy.lambdify([var], N.origin.locate_new('na', curve_diff).express_coordinates(N))
|
||||
# And a torsion function:
|
||||
torsion_fn = sympy.lambdify([var], torsion(curve_expr, var))
|
||||
# Generate values of 'var' to use:
|
||||
vs = numpy.arange(a, b, step)
|
||||
# Evaluate derivative function & torsion function over these:
|
||||
d = numpy.linalg.norm(numpy.array(diff_fn(vs)).T, axis=1)
|
||||
torsions = torsion_fn(vs)
|
||||
# Turn this into basically a left Riemann sum (I'm lazy):
|
||||
return -(d * torsions * step).sum()
|
||||
@ -1,24 +0,0 @@
|
||||
import numpy
|
||||
import quaternion
|
||||
|
||||
def conjugate_by(vec, quat):
|
||||
"""Turn 'vec' to a quaternion, conjugate it by 'quat', and return it."""
|
||||
q2 = quat * vec2quat(vec) * quat.conjugate()
|
||||
return quaternion.as_float_array(q2)[:,1:]
|
||||
|
||||
def rotation_quaternion(axis, angle):
|
||||
"""Returns a quaternion for rotating by some axis and angle.
|
||||
|
||||
Inputs:
|
||||
axis -- numpy array of shape (3,), with axis to rotate around
|
||||
angle -- angle in radians by which to rotate
|
||||
"""
|
||||
qc = numpy.cos(angle / 2)
|
||||
qs = numpy.sin(angle / 2)
|
||||
qv = qs * axis
|
||||
return numpy.quaternion(qc, qv[0], qv[1], qv[2])
|
||||
|
||||
def vec2quat(vs):
|
||||
qs = numpy.zeros(vs.shape[0], dtype=numpy.quaternion)
|
||||
quaternion.as_float_array(qs)[:,1:4] = vs
|
||||
return qs
|
||||
@ -1,232 +0,0 @@
|
||||
#!/usr/bin/env python3
|
||||
|
||||
import sys
|
||||
import numpy
|
||||
import stl.mesh
|
||||
|
||||
# TODO:
|
||||
# - This still has some strange errors around high curvature.
|
||||
# They are plainly visible in the cross-section.
|
||||
# - Check rotation direction
|
||||
# - Fix phase, which only works if 0 (due to how I work with y)
|
||||
# Things don't seem to line up right.
|
||||
# - Why is there still a gap when using Array modifier?
|
||||
# Check beginning and ending vertices maybe
|
||||
# - Organize this so that it generates both meshes when run
|
||||
|
||||
# This is all rather tightly-coupled. Almost everything is specific
|
||||
# to the isosurface I was trying to generate. walk_curve may be able
|
||||
# to generalize to some other shapes.
|
||||
class ExplicitSurfaceThing(object):
|
||||
def __init__(self, freq, phase, scale, inner, outer, rad, ext_phase):
|
||||
self.freq = freq
|
||||
self.phase = phase
|
||||
self.scale = scale
|
||||
self.inner = inner
|
||||
self.outer = outer
|
||||
self.rad = rad
|
||||
self.ext_phase = ext_phase
|
||||
|
||||
def angle(self, z):
|
||||
return self.freq*z + self.phase
|
||||
|
||||
def max_z(self):
|
||||
# This value is the largest |z| for which 'radical' >= 0
|
||||
# (thus, for x_cross to have a valid solution)
|
||||
return (numpy.arcsin(self.rad / self.inner) - self.phase) / self.freq
|
||||
|
||||
def radical(self, z):
|
||||
return self.rad*self.rad - self.inner*self.inner * (numpy.sin(self.angle(z)))**2
|
||||
|
||||
# Implicit curve function
|
||||
def F(self, x, z):
|
||||
return (self.outer*x - self.inner*numpy.cos(self.angle(z)))**2 + (self.inner*numpy.sin(self.angle(z)))**2 - self.rad**2
|
||||
|
||||
# Partial 1st derivatives of F:
|
||||
def F_x(self, x, z):
|
||||
return 2 * self.outer * self.outer * x - 2 * self.inner * self.outer * numpy.cos(self.angle(z))
|
||||
|
||||
def F_z(self, x, z):
|
||||
return 2 * self.freq * self.inner * self.outer * numpy.sin(self.angle(z))
|
||||
|
||||
# Curvature:
|
||||
def K(self, x, z):
|
||||
a1 = self.outer**2
|
||||
a2 = x**2
|
||||
a3 = self.freq*z + self.phase
|
||||
a4 = numpy.cos(a3)
|
||||
a5 = self.inner**2
|
||||
a6 = a4**2
|
||||
a7 = self.freq**2
|
||||
a8 = numpy.sin(a3)**2
|
||||
a9 = self.outer**3
|
||||
a10 = self.inner**3
|
||||
return -((2*a7*a10*self.outer*x*a4 + 2*a7*a5*a1*a2)*a8 + (2*a7*self.inner*a9*x**3 + 2*a7*a10*self.outer*x)*a4 - 4*a7*a5*a1*a2) / ((a7*a5*a2*a8 + a5*a6 - 2*self.inner*self.outer*x*a4 + a1*a2) * numpy.sqrt(4*a7*a5*a1*a2*a8 + 4*a5*a1*a6 - 8*self.inner*a9*x*a4 + 4*a2*self.outer**4))
|
||||
|
||||
def walk_curve(self, x0, z0, eps, thresh = 1e-3, gd_thresh = 1e-7):
|
||||
x, z = x0, z0
|
||||
eps2 = eps*eps
|
||||
verts = []
|
||||
iters = 0
|
||||
# Until we return to the same point at which we started...
|
||||
while True:
|
||||
iters += 1
|
||||
verts.append([x, 0, z])
|
||||
# ...walk around the curve by stepping perpendicular to the
|
||||
# gradient by 'eps'. So, first find the gradient:
|
||||
dx = self.F_x(x, z)
|
||||
dz = self.F_z(x, z)
|
||||
# Normalize it:
|
||||
f = 1/numpy.sqrt(dx*dx + dz*dz)
|
||||
nx, nz = dx*f, dz*f
|
||||
# Find curvature at this point because it tells us a little
|
||||
# about how far we can safely move:
|
||||
K_val = abs(self.K(x, z))
|
||||
eps_corr = 2 * numpy.sqrt(2*eps/K_val - eps*eps)
|
||||
# Scale by 'eps' and use (-dz, dx) as perpendicular:
|
||||
px, pz = -nz*eps_corr, nx*eps_corr
|
||||
# Walk in that direction:
|
||||
x += px
|
||||
z += pz
|
||||
# Moving in that direction is only good locally, and we may
|
||||
# have deviated off the curve slightly. The implicit function
|
||||
# tells us (sort of) how far away we are, and the gradient
|
||||
# tells us how to minimize that:
|
||||
#print("W: x={} z={} dx={} dz={} px={} pz={} K={} eps_corr={}".format(
|
||||
# x, z, dx, dz, px, pz, K_val, eps_corr))
|
||||
F_val = self.F(x, z)
|
||||
count = 0
|
||||
while abs(F_val) > gd_thresh:
|
||||
count += 1
|
||||
dx = self.F_x(x, z)
|
||||
dz = self.F_z(x, z)
|
||||
f = 1/numpy.sqrt(dx*dx + dz*dz)
|
||||
nx, nz = dx*f, dz*f
|
||||
# If F is negative, we want to increase it (thus, follow
|
||||
# gradient). If F is positive, we want to decrease it
|
||||
# (thus, opposite of gradient).
|
||||
F_val = self.F(x, z)
|
||||
x += -F_val*nx
|
||||
z += -F_val*nz
|
||||
# Yes, this is inefficient gradient-descent...
|
||||
diff = numpy.sqrt((x-x0)**2 + (z-z0)**2)
|
||||
#print("{} gradient-descent iters. diff = {}".format(count, diff))
|
||||
if iters > 100 and diff < thresh:
|
||||
#print("diff < eps, quitting")
|
||||
#verts.append([x, 0, z])
|
||||
break
|
||||
data = numpy.array(verts)
|
||||
return data
|
||||
|
||||
def x_cross(self, z, sign):
|
||||
# Single cross-section point in XZ for y=0. Set sign for positive
|
||||
# or negative solution.
|
||||
n1 = numpy.sqrt(self.radical(z))
|
||||
n2 = self.inner * numpy.cos(self.angle(z))
|
||||
if sign > 0:
|
||||
return (n2-n1) / self.outer
|
||||
else:
|
||||
return (n2+n1) / self.outer
|
||||
|
||||
def turn(self, points, dz):
|
||||
# Note one full revolution is dz = 2*pi/freq
|
||||
# How far to turn in radians (determined by dz):
|
||||
rad = self.angle(dz)
|
||||
c, s = numpy.cos(rad), numpy.sin(rad)
|
||||
mtx = numpy.array([
|
||||
[ c, s, 0],
|
||||
[-s, c, 0],
|
||||
[ 0, 0, 1],
|
||||
])
|
||||
return points.dot(mtx) + [0, 0, dz]
|
||||
|
||||
def screw_360(self, z0_period_start, x_init, z_init, eps, dz, thresh, endcaps=False):
|
||||
#z0 = -10 * 2*numpy.pi/freq / 2
|
||||
z0 = z0_period_start * 2*numpy.pi/self.freq / 2
|
||||
z1 = z0 + 2*numpy.pi/self.freq
|
||||
#z1 = 5 * 2*numpy.pi/freq / 2
|
||||
#z0 = 0
|
||||
#z1 = 2*numpy.pi/freq
|
||||
init_xsec = self.walk_curve(x_init, z_init, eps, thresh)
|
||||
num_xsec_steps = init_xsec.shape[0]
|
||||
zs = numpy.append(numpy.arange(z0, z1, dz), z1)
|
||||
num_screw_steps = len(zs)
|
||||
vecs = num_xsec_steps * num_screw_steps * 2
|
||||
offset = 0
|
||||
if endcaps:
|
||||
offset = num_xsec_steps
|
||||
vecs += 2*num_xsec_steps
|
||||
print("Generating {} vertices...".format(vecs))
|
||||
data = numpy.zeros(vecs, dtype=stl.mesh.Mesh.dtype)
|
||||
v = data["vectors"]
|
||||
# First endcap:
|
||||
if endcaps:
|
||||
center = init_xsec.mean(0)
|
||||
for i in range(num_xsec_steps):
|
||||
v[i][0,:] = init_xsec[(i + 1) % num_xsec_steps,:]
|
||||
v[i][1,:] = init_xsec[i,:]
|
||||
v[i][2,:] = center
|
||||
# Body:
|
||||
verts = init_xsec
|
||||
for i,z in enumerate(zs):
|
||||
verts_last = verts
|
||||
verts = self.turn(init_xsec, z-z0)
|
||||
if i > 0:
|
||||
for j in range(num_xsec_steps):
|
||||
# Vertex index:
|
||||
vi = offset + (i-1)*num_xsec_steps*2 + j*2
|
||||
v[vi][0,:] = verts[(j + 1) % num_xsec_steps,:]
|
||||
v[vi][1,:] = verts[j,:]
|
||||
v[vi][2,:] = verts_last[j,:]
|
||||
#print("Write vertex {}".format(vi))
|
||||
v[vi+1][0,:] = verts_last[(j + 1) % num_xsec_steps,:]
|
||||
v[vi+1][1,:] = verts[(j + 1) % num_xsec_steps,:]
|
||||
v[vi+1][2,:] = verts_last[j,:]
|
||||
#print("Write vertex {} (2nd half)".format(vi+1))
|
||||
# Second endcap:
|
||||
if endcaps:
|
||||
center = verts.mean(0)
|
||||
for i in range(num_xsec_steps):
|
||||
vi = num_xsec_steps * num_screw_steps * 2 + num_xsec_steps + i
|
||||
v[vi][0,:] = center
|
||||
v[vi][1,:] = verts[i,:]
|
||||
v[vi][2,:] = verts[(i + 1) % num_xsec_steps,:]
|
||||
v[:, :, 2] += z0 + self.ext_phase / self.freq
|
||||
v[:, :, :] /= self.scale
|
||||
mesh = stl.mesh.Mesh(data, remove_empty_areas=False)
|
||||
print("Beginning z: {}".format(z0/self.scale))
|
||||
print("Ending z: {}".format(z1/self.scale))
|
||||
print("Period: {}".format((z1-z0)/self.scale))
|
||||
return mesh
|
||||
|
||||
surf1 = ExplicitSurfaceThing(
|
||||
freq = 20,
|
||||
phase = 0,
|
||||
scale = 1/16, # from libfive
|
||||
inner = 0.4 * 1/16,
|
||||
outer = 2.0 * 1/16,
|
||||
rad = 0.3 * 1/16,
|
||||
ext_phase = 0)
|
||||
|
||||
z_init = 0
|
||||
x_init = surf1.x_cross(z_init, 1)
|
||||
mesh1 = surf1.screw_360(-10, x_init, z_init, 0.000002, 0.001, 5e-4)
|
||||
fname = "spiral_inner0_one_period.stl"
|
||||
print("Writing {}...".format(fname))
|
||||
mesh1.save(fname)
|
||||
|
||||
surf2 = ExplicitSurfaceThing(
|
||||
freq = 10,
|
||||
phase = 0,
|
||||
scale = 1/16, # from libfive
|
||||
inner = 0.9 * 1/16,
|
||||
outer = 2.0 * 1/16,
|
||||
rad = 0.3 * 1/16,
|
||||
ext_phase = numpy.pi/2)
|
||||
|
||||
z_init = 0
|
||||
x_init = surf2.x_cross(z_init, 1)
|
||||
mesh2 = surf2.screw_360(-5, x_init, z_init, 0.000002, 0.001, 5e-4)
|
||||
fname = "spiral_outer90_one_period.stl"
|
||||
print("Writing {}...".format(fname))
|
||||
mesh2.save(fname)
|
||||
@ -8,20 +8,20 @@
|
||||
- https://en.wikipedia.org/wiki/Polygon_triangulation - do this to
|
||||
fix my wave example!
|
||||
- http://www.polygontriangulation.com/2018/07/triangulation-algorithm.html
|
||||
- Clean up `examples.ram_horn_branch()`. The way I clean it up might
|
||||
- Clean up examples.ram_horn_branch(). The way I clean it up might
|
||||
help inform some cleaner designs.
|
||||
- I really need to standardize some of the behavior of fundamental
|
||||
operations (with regard to things like sizes they generate). This
|
||||
is behavior that, if it changes, will change a lot of things that I'm
|
||||
trying to keep consistent so that my examples still work.
|
||||
- Winding order. It is consistent through seemingly everything,
|
||||
except for reflection and `close_boundary_simple`. (When there are
|
||||
two parallel boundaries joined with something like
|
||||
`join_boundary_simple`, traversing these boundaries in their actual
|
||||
order to generate triangles - like in `close_boundary_simple` - will
|
||||
produce opposite winding order on each. Imagine a transparent clock:
|
||||
seen from the front, it moves clockwise, but seen from the back, it
|
||||
moves counter-clockwise.)
|
||||
- Winding order. It is consistent through seemingly
|
||||
everything, except for reflection and close_boundary_simple.
|
||||
(When there are two parallel boundaries joined with something like
|
||||
join_boundary_simple, traversing these boundaries in their actual order
|
||||
to generate triangles - like in close_boundary_simple - will produce
|
||||
opposite winding order on each. Imagine a transparent clock: seen from the
|
||||
front, it moves clockwise, but seen from the back, it moves
|
||||
counter-clockwise.)
|
||||
- File that bug that I've seen in trimesh/three.js
|
||||
(see trimesh_fail.ipynb)
|
||||
- Why do I get the weird zig-zag pattern on the triangles,
|
||||
@ -39,7 +39,7 @@
|
||||
does seem to turn 'error' just to noise, and in its own way this
|
||||
is preferable.
|
||||
- Integrate parallel_transport work and reuse what I can
|
||||
- `/mnt/dev/graphics_misc/isosurfaces_2018_2019` - perhaps include my
|
||||
- /mnt/dev/graphics_misc/isosurfaces_2018_2019 - perhaps include my
|
||||
spiral isosurface stuff from here
|
||||
|
||||
## Abstractions
|
||||
@ -81,4 +81,4 @@
|
||||
|
||||
## Other thoughts
|
||||
|
||||
- Why do I never use the term "extruding" to describe what I'm doing?
|
||||
- Why do I never use the term "extruding" to describe what I'm doing?
|
||||
@ -1,10 +0,0 @@
|
||||
#!/bin/sh
|
||||
|
||||
g++ -L/usr/local/Cellar/cgal/4.12/lib -I/usr/local/Cellar/cgal/4.12/include \
|
||||
-L/usr/local/Cellar/gmp/6.1.2_2/lib \
|
||||
-L/usr/local/Cellar/mpfr/4.0.1/lib \
|
||||
-L/usr/local/Cellar/boost/1.67.0_1/lib \
|
||||
-DCGAL_CONCURRENT_MESH_3 \
|
||||
-lCGAL -lgmp -lmpfr -lboost_thread-mt \
|
||||
./mesh_an_implicit_function.cpp \
|
||||
-o mesh_an_implicit_function.o
|
||||
@ -1,173 +0,0 @@
|
||||
#+TITLE: CGAL dabbling
|
||||
#+DATE: <2018-08-06 Mon>
|
||||
#+AUTHOR: Hodapp
|
||||
#+EMAIL: hodapp87@gmail.com
|
||||
#+OPTIONS: ':nil *:t -:t ::t <:t H:3 \n:nil ^:t arch:headline
|
||||
#+OPTIONS: author:t c:nil creator:comment d:(not "LOGBOOK") date:t
|
||||
#+OPTIONS: e:t email:nil f:t inline:t num:t p:nil pri:nil stat:t
|
||||
#+OPTIONS: tags:t tasks:t tex:t timestamp:t toc:t todo:t |:t
|
||||
#+DESCRIPTION:
|
||||
#+EXCLUDE_TAGS: noexport
|
||||
#+KEYWORDS:
|
||||
#+LANGUAGE: en
|
||||
#+SELECT_TAGS: export
|
||||
|
||||
# By default I do not want that source code blocks are evaluated on export. Usually
|
||||
# I want to evaluate them interactively and retain the original results.
|
||||
#+PROPERTY: header-args :eval never-export :export both
|
||||
|
||||
- CGAL is one of the most insanely cryptic and impenetrable libraries
|
||||
I have found.
|
||||
- Where I am stuck now:
|
||||
- I can use 1D features in
|
||||
[[https://doc.cgal.org/latest/Mesh_3/Mesh_3_2mesh_two_implicit_spheres_with_balls_8cpp-example.html][Mesh_3/mesh_two_implicit_spheres_with_balls.cpp]] but this is the
|
||||
wrong sort of data (it's a 3D mesh, yes, but not a surface mesh)
|
||||
and "medit" and "vtu" are all it can write. How do I extract a
|
||||
surface mesh to get a Polyhedron_3 so that I can write an OBJ?
|
||||
- Or: Is there any way to use 1D features with a surface mesh?
|
||||
- Even if I manage to put this Mesh_3 into the right form I have no
|
||||
guarantees that it is a good surface mesh. It uses a different
|
||||
algorithm than for the surface mesh - I'm simply trying to take
|
||||
the surface of the result. I should also expect this algorithm to
|
||||
be something more like O(N^3) with mesh resolution because it must
|
||||
fill the volume with tetrahedrons, and I then just throw away all
|
||||
of these.
|
||||
- This is based around the source for [[https://doc.cgal.org/latest/Mesh_3/Mesh_3_2mesh_implicit_sphere_8cpp-example.html][mesh_implicit_sphere.cpp]] from
|
||||
[[https://doc.cgal.org/latest/Mesh_3/index.html][CGAL: 3D Mesh Generation]].
|
||||
- Why am I not using [[https://doc.cgal.org/latest/Mesh_3/group__PkgMesh__3Functions.html#ga68ca38989087644fb6b893eb566be9ea][facets_in_complex_3_to_triangle_mesh()]]?
|
||||
- [[https://doc.cgal.org/latest/Mesh_3/Mesh_3_2mesh_two_implicit_spheres_with_balls_8cpp-example.html][Mesh_3/mesh_two_implicit_spheres_with_balls.cpp]] shows the use of
|
||||
[[https://doc.cgal.org/latest/Mesh_3/classCGAL_1_1Mesh__domain__with__polyline__features__3.html][Mesh_domain_with_polyline_features_3]] to explicitly give features to
|
||||
preserve
|
||||
- This uses the ~make_mesh_3~ call while I'm using
|
||||
~make_surface_mesh~...
|
||||
- I can only call ~CGAL::print_polyhedron_wavefront~ if I have a
|
||||
~Polyhedron~ and I cannot figure out how to get one from a ~C3t3~.
|
||||
So, perhaps I am stuck with this "medit" mesh format, or else VTU.
|
||||
- I think the terminology is trying to tell me: this isn't a surface
|
||||
mesh, it's a 3D mesh made of tetrahedra. It sounds like the type
|
||||
of data from a Delaunay triangulation is inherently rather
|
||||
different from a surface mesh.
|
||||
- https://doc.cgal.org/latest/Mesh_3/index.html#title24 does the
|
||||
opposite direction of what I need
|
||||
- There are [[https://doc.cgal.org/latest/Polygon_mesh_processing/group__PMP__detect__features__grp.html][Feature Detection Functions]], but ~detect_features~ is in
|
||||
[[https://doc.cgal.org/latest/Mesh_3/classCGAL_1_1Polyhedral__complex__mesh__domain__3.html][Polyhedral_complex_mesh_domain_3]] and
|
||||
[[https://doc.cgal.org/latest/Mesh_3/classCGAL_1_1Polyhedral__mesh__domain__with__features__3.html#a5a868ac7b8673436766d28b7a80d2826][Polyhedral_mesh_domain_with_features_3]]
|
||||
- Domains:
|
||||
- [[https://doc.cgal.org/latest/Mesh_3/classMeshDomain__3.html][MeshDomain_3]] concept
|
||||
- [[https://doc.cgal.org/latest/Mesh_3/classMeshDomainWithFeatures__3.html][MeshDomainWithFeatures_3]] concept
|
||||
- Why make_surface_mesh, Implicit_surface_3,
|
||||
Complex_2_in_triangulation_3 vs. make_mesh_3,
|
||||
Implicit_mesh_domain_3, Mesh_complex_3_in_triangulation_3?
|
||||
- [[https://doc.cgal.org/latest/Mesh_3/group__PkgMesh3Parameters.html#ga0a990b28d55157c62d4bfd2624d408af][parameters::features()]] - can do 1-dimensional features
|
||||
- now how do I use it?
|
||||
- [[https://doc.cgal.org/latest/Surface_mesher/index.html][make_surface_mesh]] args:
|
||||
- ~SurfaceMeshC2T3& c2t3~ - must be a model of concept [[https://doc.cgal.org/latest/Surface_mesher/classSurfaceMeshComplex__2InTriangulation__3.html][SurfaceMeshComplex_2InTriangulation_3]]
|
||||
- ~SurfaceMeshTraits::Surface_3 surface~ - must be a model of concept [[https://doc.cgal.org/latest/Surface_mesher/classSurface__3.html][Surface_3]]
|
||||
- ~SurfaceMeshTraits traits~ - must be a model of concept [[https://doc.cgal.org/latest/Surface_mesher/classSurfaceMeshTraits__3.html][SurfaceMeshTraits_3]]
|
||||
- ~FacetsCriteria criteria~ - must be a model of concept [[https://doc.cgal.org/latest/Surface_mesher/classSurfaceMeshFacetsCriteria__3.html][SurfaceMeshFacetsCriteria_3]]
|
||||
- ~Tag~
|
||||
- "This algorithm of CGAL::make_surface_mesh is designed for smooth
|
||||
implicit surfaces. If your implicit surface is not smooth, then the
|
||||
sharp features of the surface will not be meshed correctly."
|
||||
|
||||
#+BEGIN_SRC elisp
|
||||
(setq org-confirm-babel-evaluate nil)
|
||||
(setq org-src-fontify-natively t)
|
||||
(setq org-src-tab-acts-natively t)
|
||||
(org-version)
|
||||
#+END_SRC
|
||||
|
||||
#+RESULTS:
|
||||
: 9.1.9
|
||||
|
||||
#+BEGIN_SRC sh :results verbatim
|
||||
gcc --version
|
||||
#+END_SRC
|
||||
|
||||
#+RESULTS:
|
||||
: Apple LLVM version 9.1.0 (clang-902.0.39.2)
|
||||
: Target: x86_64-apple-darwin17.7.0
|
||||
: Thread model: posix
|
||||
: InstalledDir: /Library/Developer/CommandLineTools/usr/bin
|
||||
|
||||
#+NAME: includes
|
||||
#+BEGIN_SRC C++
|
||||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
|
||||
#include <CGAL/Mesh_triangulation_3.h>
|
||||
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
|
||||
#include <CGAL/Mesh_criteria_3.h>
|
||||
|
||||
#include <CGAL/Implicit_mesh_domain_3.h>
|
||||
#include <CGAL/make_mesh_3.h>
|
||||
#+END_SRC
|
||||
|
||||
#+NAME: typesDomain
|
||||
#+BEGIN_SRC C++
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||||
typedef K::FT FT;
|
||||
typedef K::Point_3 Point;
|
||||
typedef FT (Function)(const Point&);
|
||||
typedef CGAL::Implicit_mesh_domain_3<Function,K> Mesh_domain;
|
||||
|
||||
#ifdef CGAL_CONCURRENT_MESH_3
|
||||
typedef CGAL::Parallel_tag Concurrency_tag;
|
||||
#else
|
||||
typedef CGAL::Sequential_tag Concurrency_tag;
|
||||
#endif
|
||||
#+END_SRC
|
||||
|
||||
#+NAME: typesTriangulation
|
||||
#+BEGIN_SRC C++
|
||||
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
|
||||
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
|
||||
#+END_SRC
|
||||
|
||||
#+NAME: typesCriteria
|
||||
#+BEGIN_SRC C++
|
||||
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
|
||||
#+END_SRC
|
||||
|
||||
#+NAME: sphereFunction
|
||||
#+BEGIN_SRC C++
|
||||
FT sphere_function (const Point& p) {
|
||||
return CGAL::squared_distance(p, Point(CGAL::ORIGIN))-1;
|
||||
}
|
||||
#+END_SRC
|
||||
|
||||
#+BEGIN_SRC C++ :noweb yes :flags -L/usr/local/Cellar/cgal/4.12/lib -I/usr/local/Cellar/cgal/4.12/include -L/usr/local/Cellar/gmp/6.1.2_2/lib -L/usr/local/Cellar/mpfr/4.0.1/lib -L/usr/local/Cellar/boost/1.67.0_1/lib -DCGAL_CONCURRENT_MESH_3 -lCGAL -lgmp -lmpfr -lboost_thread-mt
|
||||
<<includes>>
|
||||
|
||||
using namespace CGAL::parameters;
|
||||
|
||||
<<typesDomain>>
|
||||
<<typesTriangulation>>
|
||||
<<typesCriteria>>
|
||||
<<sphereFunction>>
|
||||
|
||||
int main()
|
||||
{
|
||||
// Domain (Warning: Sphere_3 constructor uses squared radius !)
|
||||
Mesh_domain domain(sphere_function,
|
||||
K::Sphere_3(CGAL::ORIGIN, 2.));
|
||||
// Mesh criteria
|
||||
Mesh_criteria criteria(facet_angle=30, facet_size=0.1, facet_distance=0.025,
|
||||
cell_radius_edge_ratio=2, cell_size=0.1);
|
||||
|
||||
std::cout << "Generating..." << std::endl;
|
||||
// Mesh generation
|
||||
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
|
||||
|
||||
// Output
|
||||
std::ofstream medit_file("out.mesh");
|
||||
c3t3.output_to_medit(medit_file);
|
||||
|
||||
std::cout << "Done" << std::endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
#+END_SRC
|
||||
|
||||
#+RESULTS:
|
||||
| Generating... |
|
||||
| Done |
|
||||
Binary file not shown.
Binary file not shown.
@ -1,171 +0,0 @@
|
||||
// Taken from:
|
||||
// https://doc.cgal.org/latest/Surface_mesher/Surface_mesher_2mesh_an_implicit_function_8cpp-example.html
|
||||
// https://doc.cgal.org/latest/Surface_mesher/index.html
|
||||
// https://doc.cgal.org/latest/Mesh_3/index.html
|
||||
|
||||
#include <CGAL/Surface_mesh_default_triangulation_3.h>
|
||||
#include <CGAL/Complex_2_in_triangulation_3.h>
|
||||
#include <CGAL/IO/Complex_2_in_triangulation_3_file_writer.h>
|
||||
#ifdef CGAL_FACETS_IN_COMPLEX_2_TO_TRIANGLE_MESH_H
|
||||
#include <CGAL/IO/facets_in_complex_2_to_triangle_mesh.h>
|
||||
#else
|
||||
// NixOS currently has CGAL 4.11, not 4.12. I guess 4.12 is needed
|
||||
// for this. The #ifdef is unnecessary, but the header and call for
|
||||
// below are deprecated.
|
||||
#include <CGAL/IO/output_surface_facets_to_polyhedron.h>
|
||||
#endif
|
||||
#include <CGAL/make_surface_mesh.h>
|
||||
#include <CGAL/Surface_mesh.h>
|
||||
#include <CGAL/Implicit_surface_3.h>
|
||||
#include <CGAL/IO/print_wavefront.h>
|
||||
#include <CGAL/Polyhedron_3.h>
|
||||
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <limits>
|
||||
#include <algorithm>
|
||||
|
||||
// Triangulation
|
||||
typedef CGAL::Surface_mesh_default_triangulation_3 Tr;
|
||||
typedef CGAL::Complex_2_in_triangulation_3<Tr> C2t3;
|
||||
|
||||
// Domain?
|
||||
typedef Tr::Geom_traits GT;
|
||||
typedef GT::Sphere_3 Sphere_3;
|
||||
typedef GT::Point_3 Point_3;
|
||||
typedef GT::Vector_3 Vector_3;
|
||||
typedef GT::FT FT;
|
||||
typedef FT (*Function)(Point_3);
|
||||
typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;
|
||||
// how does this differ from CGAL::Implicit_mesh_domain_3<Function,K>?
|
||||
|
||||
typedef CGAL::Polyhedron_3<GT> Polyhedron;
|
||||
|
||||
FT sphere_function(Point_3 p) {
|
||||
Point_3 p2(p.x() + 0.1 * cos(p.x() * 20),
|
||||
p.y(),
|
||||
p.z() + 0.1 * sin(p.z() * 4));
|
||||
const FT x2=p2.x()*p2.x(), y2=p2.y()*p2.y(), z2=p2.z()*p2.z();
|
||||
return x2+y2+z2-1;
|
||||
}
|
||||
|
||||
Vector_3 sphere_gradient(Point_3 p) {
|
||||
float A = 0.1;
|
||||
float B = 0.1;
|
||||
float F1 = 20;
|
||||
float F2 = 4;
|
||||
return Vector_3(2*(A*cos(p.x()*F1) + p.x())*(1 - A*F1*sin(p.x()*F1)),
|
||||
2*p.y(),
|
||||
2*(B*sin(p.z()*F2) + p.z())*(1 + B*F2*cos(p.z()*F2)));
|
||||
}
|
||||
|
||||
FT spiral_function(Point_3 p) {
|
||||
float outer = 2.0;
|
||||
float inner = 0.4; // 0.9
|
||||
float freq = 20; // 10
|
||||
float phase = M_PI; // 3 * M_PI / 2;
|
||||
float thresh = 0.3;
|
||||
const FT d1 = p.y()*outer - inner * sin(p.x()*freq + phase);
|
||||
const FT d2 = p.z()*outer - inner * cos(p.x()*freq + phase);
|
||||
return std::max(sqrt(d1*d1 + d2*d2) - thresh,
|
||||
p.x()*p.x() + p.y()*p.y() + p.z()*p.z() - 1.9*1.9);
|
||||
}
|
||||
|
||||
Vector_3 spiral_gradient(Point_3 p) {
|
||||
float outer = 2.0;
|
||||
float inner = 0.4;
|
||||
float freq = 20;
|
||||
float phase = M_PI;
|
||||
float thresh = 0.3;
|
||||
// "block([%1,%2,%3,%4,%5,%6],%1:P+x*F,%2:cos(%1),%3:z*O-I*%2,%4:sin(%1),%5:y*O-I*%4,%6:1/sqrt(%5^2+%3^2),[((2*F*I*%3*%4-2*F*I*%2*%5)*%6)/2,O*%5*%6,O*%3*%6])"
|
||||
float v1 = phase + p.x() * freq;
|
||||
float v2 = cos(v1);
|
||||
float v3 = p.z() * outer - inner * v2;
|
||||
float v4 = sin(v1);
|
||||
float v5 = p.y() * outer - inner * v4;
|
||||
float v6 = 1.0 / sqrt(v5*v5 + v3*v3);
|
||||
return Vector_3(((2*freq*inner*v3*v4-2*freq*inner*v2*v5)*v6)/2,
|
||||
outer * v5 * v6,
|
||||
outer * v3 * v6);
|
||||
}
|
||||
|
||||
int main() {
|
||||
Tr tr; // 3D-Delaunay triangulation
|
||||
C2t3 c2t3 (tr); // 2D-complex in 3D-Delaunay triangulation
|
||||
|
||||
FT bounding_sphere_rad = 2.0;
|
||||
|
||||
// defining the surface
|
||||
Surface_3 surface(spiral_function, // pointer to function
|
||||
Sphere_3(CGAL::ORIGIN,
|
||||
bounding_sphere_rad*bounding_sphere_rad)); // bounding sphere
|
||||
|
||||
std::string fname("spiral_thing4.obj");
|
||||
float angular_bound = 30;
|
||||
float radius_bound = 0.02;
|
||||
float distance_bound = 0.02;
|
||||
|
||||
// Note that "2." above is the *squared* radius of the bounding sphere!
|
||||
// defining meshing criteria
|
||||
CGAL::Surface_mesh_default_criteria_3<Tr> criteria(
|
||||
angular_bound, radius_bound, distance_bound);
|
||||
|
||||
std::cout << "angular bound = " << angular_bound << ", "
|
||||
<< "radius bound = " << radius_bound << ", "
|
||||
<< "distance bound = " << distance_bound << std::endl;
|
||||
|
||||
std::cout << "Making surface mesh..." << std::endl;
|
||||
// meshing surface
|
||||
CGAL::make_surface_mesh(c2t3, surface, criteria, CGAL::Manifold_tag());
|
||||
std::cout << "Vertices: " << tr.number_of_vertices() << std::endl;
|
||||
|
||||
// This didn't work on some calls instead of 'poly':
|
||||
//CGAL::Surface_mesh<Point_3> sm;
|
||||
Polyhedron poly;
|
||||
std::cout << "Turning facets to triangle mesh..." << std::endl;
|
||||
|
||||
#ifdef CGAL_FACETS_IN_COMPLEX_2_TO_TRIANGLE_MESH_H
|
||||
CGAL::facets_in_complex_2_to_triangle_mesh(c2t3, poly);
|
||||
#else
|
||||
CGAL::output_surface_facets_to_polyhedron(c2t3, poly);
|
||||
#endif
|
||||
|
||||
FT err = 0.0;
|
||||
FT inf = std::numeric_limits<FT>::infinity();
|
||||
for (Polyhedron::Point_iterator it = poly.points_begin();
|
||||
it != poly.points_end();
|
||||
++it)
|
||||
{
|
||||
|
||||
FT rate = 2e-6;
|
||||
FT f0 = abs(spiral_function(*it));
|
||||
FT f;
|
||||
for (int i = 0; i < 100; ++i) {
|
||||
f = spiral_function(*it);
|
||||
Vector_3 grad(spiral_gradient(*it));
|
||||
|
||||
*it -= grad * rate * (f > 0 ? 1 : -1);
|
||||
/*
|
||||
std::cout << "Iter " << i << ": "
|
||||
<< "F(" << it->x() << "," << it->y() << "," << it->z()
|
||||
<< ")=" << f << ", F'=" << grad << std::endl;
|
||||
*/
|
||||
}
|
||||
//FT diff = (abs(f) - abs(f0)) / f0;
|
||||
/*
|
||||
std::cout << "F(" << it->x() << "," << it->y() << "," << it->z()
|
||||
<< "): " << f0 << " to " << f << std::endl;
|
||||
*/
|
||||
|
||||
err += f * f;
|
||||
}
|
||||
err = sqrt(err);
|
||||
std::cout << "RMS isosurface distance: " << err << std::endl;
|
||||
|
||||
std::cout << "Mesh vertices: " << poly.size_of_vertices() << ", "
|
||||
<< "facets: " << poly.size_of_facets() << std::endl;
|
||||
|
||||
std::cout << "Writing to " << fname << "..." << std::endl;
|
||||
std::ofstream ofs(fname);
|
||||
CGAL::print_polyhedron_wavefront(ofs, poly);
|
||||
}
|
||||
@ -1,55 +0,0 @@
|
||||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
|
||||
#include <CGAL/Mesh_triangulation_3.h>
|
||||
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
|
||||
#include <CGAL/Mesh_criteria_3.h>
|
||||
|
||||
#include <CGAL/Implicit_mesh_domain_3.h>
|
||||
#include <CGAL/make_mesh_3.h>
|
||||
|
||||
// Domain
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||||
typedef K::FT FT;
|
||||
typedef K::Point_3 Point;
|
||||
typedef FT (Function)(const Point&);
|
||||
typedef CGAL::Implicit_mesh_domain_3<Function,K> Mesh_domain;
|
||||
|
||||
#ifdef CGAL_CONCURRENT_MESH_3
|
||||
typedef CGAL::Parallel_tag Concurrency_tag;
|
||||
#else
|
||||
typedef CGAL::Sequential_tag Concurrency_tag;
|
||||
#endif
|
||||
|
||||
// Triangulation
|
||||
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
|
||||
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
|
||||
|
||||
// Criteria
|
||||
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
|
||||
|
||||
// To avoid verbose function and named parameters call
|
||||
using namespace CGAL::parameters;
|
||||
|
||||
// Function
|
||||
FT sphere_function (const Point& p)
|
||||
{ return CGAL::squared_distance(p, Point(CGAL::ORIGIN))-1; }
|
||||
|
||||
int main()
|
||||
{
|
||||
// Domain (Warning: Sphere_3 constructor uses squared radius !)
|
||||
Mesh_domain domain(sphere_function,
|
||||
K::Sphere_3(CGAL::ORIGIN, 2.));
|
||||
// Mesh criteria
|
||||
Mesh_criteria criteria(facet_angle=30, facet_size=0.1, facet_distance=0.025,
|
||||
cell_radius_edge_ratio=2, cell_size=0.1);
|
||||
|
||||
// Mesh generation
|
||||
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
|
||||
|
||||
// Output
|
||||
std::ofstream medit_file("out.mesh");
|
||||
c3t3.output_to_medit(medit_file);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
@ -1,106 +0,0 @@
|
||||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
|
||||
#include <CGAL/Mesh_triangulation_3.h>
|
||||
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
|
||||
#include <CGAL/Mesh_criteria_3.h>
|
||||
|
||||
#include <CGAL/Labeled_mesh_domain_3.h>
|
||||
#include <CGAL/Mesh_domain_with_polyline_features_3.h>
|
||||
#include <CGAL/make_mesh_3.h>
|
||||
|
||||
#include <CGAL/IO/print_wavefront.h>
|
||||
|
||||
// Kernel
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||||
|
||||
// Domain
|
||||
typedef K::FT FT;
|
||||
typedef K::Point_3 Point;
|
||||
typedef FT (Function)(const Point&);
|
||||
typedef CGAL::Mesh_domain_with_polyline_features_3<
|
||||
CGAL::Labeled_mesh_domain_3<K> > Mesh_domain;
|
||||
|
||||
// Polyline
|
||||
typedef std::vector<Point> Polyline_3;
|
||||
typedef std::list<Polyline_3> Polylines;
|
||||
|
||||
#ifdef CGAL_CONCURRENT_MESH_3
|
||||
typedef CGAL::Parallel_tag Concurrency_tag;
|
||||
#else
|
||||
typedef CGAL::Sequential_tag Concurrency_tag;
|
||||
#endif
|
||||
|
||||
// Triangulation
|
||||
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
|
||||
|
||||
typedef CGAL::Mesh_complex_3_in_triangulation_3<
|
||||
Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_index> C3t3;
|
||||
|
||||
// Criteria
|
||||
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
|
||||
|
||||
// To avoid verbose function and named parameters call
|
||||
using namespace CGAL::parameters;
|
||||
|
||||
// Function
|
||||
FT sphere_function1 (const Point& p)
|
||||
{ return CGAL::squared_distance(p, Point(CGAL::ORIGIN))-2; }
|
||||
|
||||
FT sphere_function2 (const Point& p)
|
||||
{ return CGAL::squared_distance(p, Point(2, 0, 0))-2; }
|
||||
|
||||
FT sphere_function (const Point& p)
|
||||
{
|
||||
if(sphere_function1(p) < 0 || sphere_function2(p) < 0)
|
||||
return -1;
|
||||
else
|
||||
return 1;
|
||||
}
|
||||
|
||||
#include <cmath>
|
||||
|
||||
int main()
|
||||
{
|
||||
// Domain (Warning: Sphere_3 constructor uses squared radius !)
|
||||
Mesh_domain domain =
|
||||
Mesh_domain::create_implicit_mesh_domain(sphere_function,
|
||||
K::Sphere_3(Point(1, 0, 0), 6.));
|
||||
|
||||
// Mesh criteria
|
||||
Mesh_criteria criteria(edge_size = 0.15,
|
||||
facet_angle = 25, facet_size = 0.15,
|
||||
cell_radius_edge_ratio = 2, cell_size = 0.15);
|
||||
|
||||
// Create edge that we want to preserve
|
||||
Polylines polylines (1);
|
||||
Polyline_3& polyline = polylines.front();
|
||||
|
||||
for(int i = 0; i < 360; ++i)
|
||||
{
|
||||
Point p (1, std::cos(i*CGAL_PI/180), std::sin(i*CGAL_PI/180));
|
||||
polyline.push_back(p);
|
||||
}
|
||||
polyline.push_back(polyline.front()); // close the line
|
||||
|
||||
// Insert edge in domain
|
||||
domain.add_features(polylines.begin(), polylines.end());
|
||||
|
||||
// Mesh generation without feature preservation
|
||||
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
|
||||
CGAL::parameters::no_features());
|
||||
|
||||
std::ofstream medit_file("out-no-protection.vtu");
|
||||
c3t3.output_to_medit(medit_file);
|
||||
medit_file.close();
|
||||
c3t3.clear();
|
||||
|
||||
// Mesh generation with feature preservation
|
||||
c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
|
||||
|
||||
// Output
|
||||
medit_file.open("out-with-protection.vtu");
|
||||
c3t3.output_to_medit(medit_file);
|
||||
medit_file.close();
|
||||
|
||||
return 0;
|
||||
}
|
||||
@ -1,3 +0,0 @@
|
||||
#!/bin/sh
|
||||
|
||||
nix-shell -p python37Packages.numpy-stl
|
||||
@ -1,21 +0,0 @@
|
||||
#!/usr/bin/env python
|
||||
|
||||
import math
|
||||
|
||||
def f1(T, I, O, P, F):
|
||||
return lambda x: (x, (T + I * math.sin(P + x*F)) / O, (I * math.cos(P + x*F)) / O)
|
||||
|
||||
def f2(T, I, O, P, F):
|
||||
return lambda x: (x, -(T - I * math.sin(P + x*F)) / O, (I * math.cos(P + x*F)) / O)
|
||||
|
||||
f = f2(O=2.0, I=0.4, F=20, P=math.pi, T=0.3)
|
||||
print("ply")
|
||||
print("format ascii 1.0")
|
||||
r = range(-400, 400)
|
||||
print("element vertex %d" % (len(r)))
|
||||
print("property float32 x")
|
||||
print("property float32 y")
|
||||
print("property float32 z")
|
||||
print("end_header")
|
||||
for xi in r:
|
||||
print("%f %f %f" % f(float(xi) / 200))
|
||||
@ -1,9 +0,0 @@
|
||||
{ pkgs ? import <nixpkgs> {} }:
|
||||
|
||||
let python_with_deps = pkgs.python3.withPackages
|
||||
(ps: [ps.numpy ps.numpy-stl]);
|
||||
in pkgs.stdenv.mkDerivation rec {
|
||||
name = "gfx_scratch";
|
||||
|
||||
buildInputs = with pkgs; [ python_with_deps ];
|
||||
}
|
||||
@ -1,10 +0,0 @@
|
||||
{ pkgs ? import <nixpkgs> {} }:
|
||||
|
||||
let stdenv = pkgs.stdenv;
|
||||
in stdenv.mkDerivation rec {
|
||||
name = "cgal_scratch";
|
||||
|
||||
buildInputs = with pkgs; [ cgal boost gmp mpfr ];
|
||||
}
|
||||
|
||||
# g++ -lCGAL -lmpfr -lgmp mesh_an_implicit_function.cpp -o mesh_an_implicit_function.o
|
||||
@ -1,109 +0,0 @@
|
||||
#!/usr/bin/env python3
|
||||
|
||||
import sys
|
||||
import numpy
|
||||
import stl.mesh
|
||||
|
||||
# TODO:
|
||||
# - This is a very naive triangulation strategy. It needs fixing - the
|
||||
# way it handles 'flatter' areas isn't optimal at all, even if the
|
||||
# sharper areas are much better than from CGAL or libfive.
|
||||
# - Generate just part of the mesh and then copy. It is rotationally
|
||||
# symmetric, as well as translationally symmetric at its period.
|
||||
|
||||
fname = "spiral_outer0.stl"
|
||||
freq = 20
|
||||
phase = 0
|
||||
scale = 1/16 # from libfive
|
||||
inner = 0.4 * scale
|
||||
outer = 2.0 * scale
|
||||
rad = 0.3 * scale
|
||||
|
||||
angle = lambda z: freq*z + phase
|
||||
|
||||
# z starting & ending point:
|
||||
z0 = -20*scale
|
||||
z1 = 20*scale
|
||||
# Number of z divisions:
|
||||
m = 1600
|
||||
# Number of circle points:
|
||||
n = 1000
|
||||
|
||||
dz = (z1 - z0) / (m-1)
|
||||
|
||||
data = numpy.zeros((m-1)*n*2 + 2*n, dtype=stl.mesh.Mesh.dtype)
|
||||
# Vertex count:
|
||||
# From z0 to z0+dz is n circle points joined with 2 triangles to next -> n*2
|
||||
# z0+dz to z0+dz*2 is likewise... up through (m-1) of these -> (m-1)*n*2
|
||||
# Two endcaps each have circle points & center point -> 2*n
|
||||
# Thus: (m-1)*n*2 + 2*n
|
||||
v = data["vectors"]
|
||||
print("Vertex count: {}".format(m*n*2 + 2*n))
|
||||
|
||||
verts = numpy.zeros((n, 3), dtype=numpy.float32)
|
||||
|
||||
# For every z cross-section...
|
||||
for z_idx in range(m):
|
||||
#sys.stdout.write(".")
|
||||
# z value:
|
||||
z = z0 + dz*z_idx
|
||||
# Angle of center point of circle (radians):
|
||||
# (we don't actually need to normalize this)
|
||||
rad = angle(z)
|
||||
c,s = numpy.cos(rad), numpy.sin(rad)
|
||||
# Center point of circle:
|
||||
cx, cy = (inner + outer)*numpy.cos(rad), (inner + outer)*numpy.sin(rad)
|
||||
# For every division of the circular cross-section...
|
||||
if z_idx == 0:
|
||||
# Begin with z0 endcap as a special case:
|
||||
verts_last = numpy.zeros((n, 3), dtype=numpy.float32)
|
||||
verts_last[:, 0] = cx
|
||||
verts_last[:, 1] = cy
|
||||
verts_last[:, 2] = z
|
||||
else:
|
||||
verts_last = verts
|
||||
verts = numpy.zeros((n, 3), dtype=numpy.float32)
|
||||
for ang_idx in range(n):
|
||||
# Step around starting angle (the 'far' intersection of the
|
||||
# line at angle 'rad' and this circle):
|
||||
rad2 = rad + 2*ang_idx*numpy.pi/n
|
||||
# ...and generate points on the circle:
|
||||
xi = cx + outer*numpy.cos(rad2)
|
||||
yi = cy + outer*numpy.sin(rad2)
|
||||
verts[ang_idx, :] = [xi, yi, z]
|
||||
#print("i={}, z={}, rad={}, cx={}, cy={}, rad2={}, xi={}, yi={}".format(i,z,rad,cx,cy, rad2, xi, yi))
|
||||
if z_idx == 0:
|
||||
for i in range(n):
|
||||
v[i][0,:] = verts[(i + 1) % n,:]
|
||||
v[i][1,:] = verts[i,:]
|
||||
v[i][2,:] = verts_last[i,:]
|
||||
#print("Write vertex {}".format(i))
|
||||
else:
|
||||
for i in range(n):
|
||||
# Vertex index:
|
||||
vi = z_idx*n*2 + i*2 - n
|
||||
v[vi][0,:] = verts[(i + 1) % n,:]
|
||||
v[vi][1,:] = verts[i,:]
|
||||
v[vi][2,:] = verts_last[i,:]
|
||||
#print("Write vertex {}".format(vi))
|
||||
v[vi+1][0,:] = verts_last[(i + 1) % n,:]
|
||||
v[vi+1][1,:] = verts[(i + 1) % n,:]
|
||||
v[vi+1][2,:] = verts_last[i,:]
|
||||
#print("Write vertex {} (2nd half)".format(vi+1))
|
||||
|
||||
# then handle z1 endcap:
|
||||
for i in range(n):
|
||||
# See vi definition above. z_idx ends at m-1, i ends at n-1, and
|
||||
# so evaluate vi+1 (final index it wrote), add 1 for the next, and
|
||||
# then use 'i' to step one at a time:
|
||||
vi = (m-1)*n*2 + (n-1)*2 - n + 2 + i
|
||||
v[vi][0,:] = verts[i,:]
|
||||
v[vi][1,:] = verts[(i + 1) % n,:]
|
||||
v[vi][2,:] = [cx, cy, z]
|
||||
# Note winding order (1 & 2 flipped from other endcap)
|
||||
#print("Write vertex {} (endcap)".format(vi))
|
||||
|
||||
print("Writing {}...".format(fname))
|
||||
mesh = stl.mesh.Mesh(data, remove_empty_areas=False)
|
||||
mesh.save(fname)
|
||||
print("Done.")
|
||||
@ -1,201 +0,0 @@
|
||||
#!/usr/bin/env python3
|
||||
|
||||
import sys
|
||||
import numpy
|
||||
import stl.mesh
|
||||
|
||||
# TODO:
|
||||
# - Fix correction around high curvatures. It has boundary issues
|
||||
# between the functions.
|
||||
# - Check every vertex point against the actual isosurface.
|
||||
# - Check rotation direction
|
||||
# - Fix phase, which only works if 0!
|
||||
|
||||
fname = "spiral_inner0_one_period.stl"
|
||||
freq = 20
|
||||
phase = 0
|
||||
scale = 1/16 # from libfive
|
||||
inner = 0.4 * scale
|
||||
outer = 2.0 * scale
|
||||
rad = 0.3 * scale
|
||||
ext_phase = 0
|
||||
|
||||
"""
|
||||
fname = "spiral_outer90_one_period.stl"
|
||||
freq = 10
|
||||
#phase = numpy.pi/2
|
||||
phase = 0
|
||||
scale = 1/16 # from libfive
|
||||
inner = 0.9 * scale
|
||||
outer = 2.0 * scale
|
||||
rad = 0.3 * scale
|
||||
ext_phase = numpy.pi/2
|
||||
"""
|
||||
|
||||
def angle(z):
|
||||
return freq*z + phase
|
||||
|
||||
def max_z():
|
||||
# This value is the largest |z| for which 'radical' >= 0
|
||||
# (thus, for x_cross to have a valid solution)
|
||||
return (numpy.arcsin(rad / inner) - phase) / freq
|
||||
|
||||
def radical(z):
|
||||
return rad*rad - inner*inner * (numpy.sin(angle(z)))**2
|
||||
|
||||
def x_cross(z, sign):
|
||||
# Single cross-section point in XZ for y=0. Set sign for positive
|
||||
# or negative solution.
|
||||
n1 = numpy.sqrt(radical(z))
|
||||
n2 = inner * numpy.cos(angle(z))
|
||||
if sign > 0:
|
||||
return (n2-n1) / outer
|
||||
else:
|
||||
return (n2+n1) / outer
|
||||
|
||||
def curvature_cross(z, sign):
|
||||
# Curvature at a given cross-section point. This is fugly because
|
||||
# it was produced from Maxima's optimized expression.
|
||||
a1 = 1/outer
|
||||
a2 = freq**2
|
||||
a3 = phase + z*freq
|
||||
a4 = numpy.cos(a3)
|
||||
a5 = a4**2
|
||||
a6 = numpy.sin(a3)
|
||||
a7 = a6**2
|
||||
a8 = inner**2
|
||||
a9 = numpy.sqrt(rad**2 - a8*a7)
|
||||
a10 = -a2*(inner**4)*a5*a7 / (a9**3)
|
||||
a11 = 1 / a9
|
||||
a12 = -a2*a8*a5*a11
|
||||
a13 = a2*a8*a7*a11
|
||||
a14 = 1/(outer**2)
|
||||
a15 = -freq*a8*a4*a6*a11
|
||||
if sign > 0:
|
||||
return -a1*(a13+a12+a10+a2*inner*a4) / ((a14*(a15+freq*inner*a6)**2 + 1)**(3/2))
|
||||
else:
|
||||
return a1*(a13+a13+a10-a2*inner*a4) / ((a14*(a15-freq*inner*a6)**2 + 1)**(3/2))
|
||||
|
||||
def cross_section_xz(eps):
|
||||
# Generate points for a cross-section in XZ. 'eps' is the maximum
|
||||
# distance in either axis.
|
||||
verts = []
|
||||
signs = [-1, 1]
|
||||
z_start = [0, max_z()]
|
||||
z_end = [max_z(), 0]
|
||||
# Yes, this is clunky and numerical:
|
||||
for sign, z0, z1 in zip(signs, z_start, z_end):
|
||||
print("sign={} z0={} z1={}".format(sign, z0, z1))
|
||||
z = z0
|
||||
x = x_cross(z, sign)
|
||||
while (sign*z) >= (sign*z1):
|
||||
verts.append([x, 0, z])
|
||||
x_last = x
|
||||
dz = -sign*min(eps, abs(z - z1))
|
||||
if abs(dz) < 1e-8:
|
||||
break
|
||||
x = x_cross(z + dz, sign)
|
||||
#curvature = max(abs(curvature_cross(z, sign)), abs(curvature_cross(z + dz, sign)))
|
||||
curvature = abs(curvature_cross((z + dz)/2, sign))
|
||||
dx = (x - x_last) * curvature
|
||||
print("start x={} dx={} z={} dz={} curvature={}".format(x, dx, z, dz, curvature))
|
||||
while abs(dx) > eps:
|
||||
dz *= 0.8
|
||||
x = x_cross(z + dz, sign)
|
||||
curvature = abs(curvature_cross((z + dz)/2, sign))
|
||||
#curvature = max(abs(curvature_cross(z, sign)), abs(curvature_cross(z + dz, sign)))
|
||||
dx = (x - x_last) * curvature
|
||||
print("iter x={} dx={} z={} dz={} curvature={}".format(x, dx, z, dz, curvature))
|
||||
z = z + dz
|
||||
print("finish x={} z={} curvature={}".format(x, z, curvature))
|
||||
n = len(verts)
|
||||
data = numpy.zeros((n*2, 3))
|
||||
data[:n, :] = verts
|
||||
data[n:, :] = verts[::-1]
|
||||
data[n:, 2] = -data[n:, 2]
|
||||
return data
|
||||
|
||||
def turn(points, dz):
|
||||
# Note one full revolution is dz = 2*pi/freq
|
||||
# How far to turn in radians (determined by dz):
|
||||
rad = angle(dz)
|
||||
c, s = numpy.cos(rad), numpy.sin(rad)
|
||||
mtx = numpy.array([
|
||||
[ c, s, 0],
|
||||
[-s, c, 0],
|
||||
[ 0, 0, 1],
|
||||
])
|
||||
return points.dot(mtx) + [0, 0, dz]
|
||||
|
||||
def screw_360(eps, dz):
|
||||
#z0 = -10 * 2*numpy.pi/freq / 2
|
||||
z0 = -5 * 2*numpy.pi/freq / 2
|
||||
z1 = z0 + 2*numpy.pi/freq
|
||||
#z1 = 5 * 2*numpy.pi/freq / 2
|
||||
#z0 = 0
|
||||
#z1 = 2*numpy.pi/freq
|
||||
init_xsec = cross_section_xz(eps)
|
||||
num_xsec_steps = init_xsec.shape[0]
|
||||
zs = numpy.arange(z0, z1, dz)
|
||||
num_screw_steps = len(zs)
|
||||
vecs = num_xsec_steps * num_screw_steps * 2 + 2*num_xsec_steps
|
||||
print("Generating {} vertices...".format(vecs))
|
||||
data = numpy.zeros(vecs, dtype=stl.mesh.Mesh.dtype)
|
||||
v = data["vectors"]
|
||||
# First endcap:
|
||||
center = init_xsec.mean(0)
|
||||
for i in range(num_xsec_steps):
|
||||
v[i][0,:] = init_xsec[(i + 1) % num_xsec_steps,:]
|
||||
v[i][1,:] = init_xsec[i,:]
|
||||
v[i][2,:] = center
|
||||
# Body:
|
||||
verts = init_xsec
|
||||
for i,z in enumerate(zs):
|
||||
verts_last = verts
|
||||
verts = turn(init_xsec, z-z0)
|
||||
if i > 0:
|
||||
for j in range(num_xsec_steps):
|
||||
# Vertex index:
|
||||
vi = num_xsec_steps + (i-1)*num_xsec_steps*2 + j*2
|
||||
v[vi][0,:] = verts[(j + 1) % num_xsec_steps,:]
|
||||
v[vi][1,:] = verts[j,:]
|
||||
v[vi][2,:] = verts_last[j,:]
|
||||
#print("Write vertex {}".format(vi))
|
||||
v[vi+1][0,:] = verts_last[(j + 1) % num_xsec_steps,:]
|
||||
v[vi+1][1,:] = verts[(j + 1) % num_xsec_steps,:]
|
||||
v[vi+1][2,:] = verts_last[j,:]
|
||||
#print("Write vertex {} (2nd half)".format(vi+1))
|
||||
# Second endcap:
|
||||
center = verts.mean(0)
|
||||
for i in range(num_xsec_steps):
|
||||
vi = num_xsec_steps * num_screw_steps * 2 + num_xsec_steps + i
|
||||
v[vi][0,:] = center
|
||||
v[vi][1,:] = verts[i,:]
|
||||
v[vi][2,:] = verts[(i + 1) % num_xsec_steps,:]
|
||||
v[:, :, 2] += z0 + ext_phase / freq
|
||||
v[:, :, :] /= scale
|
||||
mesh = stl.mesh.Mesh(data, remove_empty_areas=False)
|
||||
print("Beginning z: {}".format(z0/scale))
|
||||
print("Ending z: {}".format(z1/scale))
|
||||
print("Period: {}".format((z1-z0)/scale))
|
||||
return mesh
|
||||
|
||||
#print("Writing {}...".format(fname))
|
||||
#mesh = stl.mesh.Mesh(data, remove_empty_areas=False)
|
||||
#mesh.save(fname)
|
||||
#print("Done.")
|
||||
|
||||
# What's up with this? Note the jump from z=0.0424 to z=0.037.
|
||||
"""
|
||||
finish x=0.13228756555322954 z=0.042403103949074046 curvature=2.451108140319032
|
||||
sign=1 z0=0.042403103949074046 z1=0
|
||||
__main__:75: RuntimeWarning: invalid value encountered in double_scalars
|
||||
start x=0.0834189730812818 dx=nan z=0.042403103949074046 dz=-0.005 curvature=nan
|
||||
finish x=0.0834189730812818 z=0.03740310394907405 curvature=nan
|
||||
"""
|
||||
# Is it because curvature is undefined there - thus the starting step
|
||||
# size of 0.005 is fine?
|
||||
|
||||
m = screw_360(0.01, 0.001)
|
||||
print("Writing {}...".format(fname))
|
||||
m.save(fname)
|
||||
@ -1,234 +0,0 @@
|
||||
#!/usr/bin/env python3
|
||||
|
||||
import sys
|
||||
import numpy
|
||||
import stl.mesh
|
||||
|
||||
# TODO:
|
||||
# - This still has some strange errors around high curvature.
|
||||
# They are plainly visible in the cross-section.
|
||||
# (errr... which errors were those? I can't see in the render)
|
||||
# - Check rotation direction
|
||||
# - Fix phase, which only works if 0 (due to how I work with y)
|
||||
# Things don't seem to line up right.
|
||||
# - Why is there still a gap when using Array modifier?
|
||||
# Check beginning and ending vertices maybe
|
||||
# - Organize this so that it generates both meshes when run
|
||||
# - Use SymPy instead of all this hard-coded stuff?
|
||||
|
||||
# This is all rather tightly-coupled. Almost everything is specific
|
||||
# to the isosurface I was trying to generate. walk_curve may be able
|
||||
# to generalize to some other shapes.
|
||||
class ExplicitSurfaceThing(object):
|
||||
def __init__(self, freq, phase, scale, inner, outer, rad, ext_phase):
|
||||
self.freq = freq
|
||||
self.phase = phase
|
||||
self.scale = scale
|
||||
self.inner = inner
|
||||
self.outer = outer
|
||||
self.rad = rad
|
||||
self.ext_phase = ext_phase
|
||||
|
||||
def angle(self, z):
|
||||
return self.freq*z + self.phase
|
||||
|
||||
def max_z(self):
|
||||
# This value is the largest |z| for which 'radical' >= 0
|
||||
# (thus, for x_cross to have a valid solution)
|
||||
return (numpy.arcsin(self.rad / self.inner) - self.phase) / self.freq
|
||||
|
||||
def radical(self, z):
|
||||
return self.rad*self.rad - self.inner*self.inner * (numpy.sin(self.angle(z)))**2
|
||||
|
||||
# Implicit curve function
|
||||
def F(self, x, z):
|
||||
return (self.outer*x - self.inner*numpy.cos(self.angle(z)))**2 + (self.inner*numpy.sin(self.angle(z)))**2 - self.rad**2
|
||||
|
||||
# Partial 1st derivatives of F:
|
||||
def F_x(self, x, z):
|
||||
return 2 * self.outer * self.outer * x - 2 * self.inner * self.outer * numpy.cos(self.angle(z))
|
||||
|
||||
def F_z(self, x, z):
|
||||
return 2 * self.freq * self.inner * self.outer * numpy.sin(self.angle(z))
|
||||
|
||||
# Curvature:
|
||||
def K(self, x, z):
|
||||
a1 = self.outer**2
|
||||
a2 = x**2
|
||||
a3 = self.freq*z + self.phase
|
||||
a4 = numpy.cos(a3)
|
||||
a5 = self.inner**2
|
||||
a6 = a4**2
|
||||
a7 = self.freq**2
|
||||
a8 = numpy.sin(a3)**2
|
||||
a9 = self.outer**3
|
||||
a10 = self.inner**3
|
||||
return -((2*a7*a10*self.outer*x*a4 + 2*a7*a5*a1*a2)*a8 + (2*a7*self.inner*a9*x**3 + 2*a7*a10*self.outer*x)*a4 - 4*a7*a5*a1*a2) / ((a7*a5*a2*a8 + a5*a6 - 2*self.inner*self.outer*x*a4 + a1*a2) * numpy.sqrt(4*a7*a5*a1*a2*a8 + 4*a5*a1*a6 - 8*self.inner*a9*x*a4 + 4*a2*self.outer**4))
|
||||
|
||||
def walk_curve(self, x0, z0, eps, thresh = 1e-3, gd_thresh = 1e-7):
|
||||
x, z = x0, z0
|
||||
eps2 = eps*eps
|
||||
verts = []
|
||||
iters = 0
|
||||
# Until we return to the same point at which we started...
|
||||
while True:
|
||||
iters += 1
|
||||
verts.append([x, 0, z])
|
||||
# ...walk around the curve by stepping perpendicular to the
|
||||
# gradient by 'eps'. So, first find the gradient:
|
||||
dx = self.F_x(x, z)
|
||||
dz = self.F_z(x, z)
|
||||
# Normalize it:
|
||||
f = 1/numpy.sqrt(dx*dx + dz*dz)
|
||||
nx, nz = dx*f, dz*f
|
||||
# Find curvature at this point because it tells us a little
|
||||
# about how far we can safely move:
|
||||
K_val = abs(self.K(x, z))
|
||||
eps_corr = 2 * numpy.sqrt(2*eps/K_val - eps*eps)
|
||||
# Scale by 'eps' and use (-dz, dx) as perpendicular:
|
||||
px, pz = -nz*eps_corr, nx*eps_corr
|
||||
# Walk in that direction:
|
||||
x += px
|
||||
z += pz
|
||||
# Moving in that direction is only good locally, and we may
|
||||
# have deviated off the curve slightly. The implicit function
|
||||
# tells us (sort of) how far away we are, and the gradient
|
||||
# tells us how to minimize that:
|
||||
#print("W: x={} z={} dx={} dz={} px={} pz={} K={} eps_corr={}".format(
|
||||
# x, z, dx, dz, px, pz, K_val, eps_corr))
|
||||
F_val = self.F(x, z)
|
||||
count = 0
|
||||
while abs(F_val) > gd_thresh:
|
||||
count += 1
|
||||
dx = self.F_x(x, z)
|
||||
dz = self.F_z(x, z)
|
||||
f = 1/numpy.sqrt(dx*dx + dz*dz)
|
||||
nx, nz = dx*f, dz*f
|
||||
# If F is negative, we want to increase it (thus, follow
|
||||
# gradient). If F is positive, we want to decrease it
|
||||
# (thus, opposite of gradient).
|
||||
F_val = self.F(x, z)
|
||||
x += -F_val*nx
|
||||
z += -F_val*nz
|
||||
# Yes, this is inefficient gradient-descent...
|
||||
diff = numpy.sqrt((x-x0)**2 + (z-z0)**2)
|
||||
#print("{} gradient-descent iters. diff = {}".format(count, diff))
|
||||
if iters > 100 and diff < thresh:
|
||||
#print("diff < eps, quitting")
|
||||
#verts.append([x, 0, z])
|
||||
break
|
||||
data = numpy.array(verts)
|
||||
return data
|
||||
|
||||
def x_cross(self, z, sign):
|
||||
# Single cross-section point in XZ for y=0. Set sign for positive
|
||||
# or negative solution.
|
||||
n1 = numpy.sqrt(self.radical(z))
|
||||
n2 = self.inner * numpy.cos(self.angle(z))
|
||||
if sign > 0:
|
||||
return (n2-n1) / self.outer
|
||||
else:
|
||||
return (n2+n1) / self.outer
|
||||
|
||||
def turn(self, points, dz):
|
||||
# Note one full revolution is dz = 2*pi/freq
|
||||
# How far to turn in radians (determined by dz):
|
||||
rad = self.angle(dz)
|
||||
c, s = numpy.cos(rad), numpy.sin(rad)
|
||||
mtx = numpy.array([
|
||||
[ c, s, 0],
|
||||
[-s, c, 0],
|
||||
[ 0, 0, 1],
|
||||
])
|
||||
return points.dot(mtx) + [0, 0, dz]
|
||||
|
||||
def screw_360(self, z0_period_start, x_init, z_init, eps, dz, thresh, endcaps=False):
|
||||
#z0 = -10 * 2*numpy.pi/freq / 2
|
||||
z0 = z0_period_start * 2*numpy.pi/self.freq / 2
|
||||
z1 = z0 + 2*numpy.pi/self.freq
|
||||
#z1 = 5 * 2*numpy.pi/freq / 2
|
||||
#z0 = 0
|
||||
#z1 = 2*numpy.pi/freq
|
||||
init_xsec = self.walk_curve(x_init, z_init, eps, thresh)
|
||||
num_xsec_steps = init_xsec.shape[0]
|
||||
zs = numpy.append(numpy.arange(z0, z1, dz), z1)
|
||||
num_screw_steps = len(zs)
|
||||
vecs = num_xsec_steps * num_screw_steps * 2
|
||||
offset = 0
|
||||
if endcaps:
|
||||
offset = num_xsec_steps
|
||||
vecs += 2*num_xsec_steps
|
||||
print("Generating {} vertices...".format(vecs))
|
||||
data = numpy.zeros(vecs, dtype=stl.mesh.Mesh.dtype)
|
||||
v = data["vectors"]
|
||||
# First endcap:
|
||||
if endcaps:
|
||||
center = init_xsec.mean(0)
|
||||
for i in range(num_xsec_steps):
|
||||
v[i][0,:] = init_xsec[(i + 1) % num_xsec_steps,:]
|
||||
v[i][1,:] = init_xsec[i,:]
|
||||
v[i][2,:] = center
|
||||
# Body:
|
||||
verts = init_xsec
|
||||
for i,z in enumerate(zs):
|
||||
verts_last = verts
|
||||
verts = self.turn(init_xsec, z-z0)
|
||||
if i > 0:
|
||||
for j in range(num_xsec_steps):
|
||||
# Vertex index:
|
||||
vi = offset + (i-1)*num_xsec_steps*2 + j*2
|
||||
v[vi][0,:] = verts[(j + 1) % num_xsec_steps,:]
|
||||
v[vi][1,:] = verts[j,:]
|
||||
v[vi][2,:] = verts_last[j,:]
|
||||
#print("Write vertex {}".format(vi))
|
||||
v[vi+1][0,:] = verts_last[(j + 1) % num_xsec_steps,:]
|
||||
v[vi+1][1,:] = verts[(j + 1) % num_xsec_steps,:]
|
||||
v[vi+1][2,:] = verts_last[j,:]
|
||||
#print("Write vertex {} (2nd half)".format(vi+1))
|
||||
# Second endcap:
|
||||
if endcaps:
|
||||
center = verts.mean(0)
|
||||
for i in range(num_xsec_steps):
|
||||
vi = num_xsec_steps * num_screw_steps * 2 + num_xsec_steps + i
|
||||
v[vi][0,:] = center
|
||||
v[vi][1,:] = verts[i,:]
|
||||
v[vi][2,:] = verts[(i + 1) % num_xsec_steps,:]
|
||||
v[:, :, 2] += z0 + self.ext_phase / self.freq
|
||||
v[:, :, :] /= self.scale
|
||||
mesh = stl.mesh.Mesh(data, remove_empty_areas=False)
|
||||
print("Beginning z: {}".format(z0/self.scale))
|
||||
print("Ending z: {}".format(z1/self.scale))
|
||||
print("Period: {}".format((z1-z0)/self.scale))
|
||||
return mesh
|
||||
|
||||
surf1 = ExplicitSurfaceThing(
|
||||
freq = 20,
|
||||
phase = 0,
|
||||
scale = 1/16, # from libfive
|
||||
inner = 0.4 * 1/16,
|
||||
outer = 2.0 * 1/16,
|
||||
rad = 0.3 * 1/16,
|
||||
ext_phase = 0)
|
||||
|
||||
z_init = 0
|
||||
x_init = surf1.x_cross(z_init, 1)
|
||||
mesh1 = surf1.screw_360(-10, x_init, z_init, 0.000002, 0.001, 5e-4)
|
||||
fname = "spiral_inner0_one_period.stl"
|
||||
print("Writing {}...".format(fname))
|
||||
mesh1.save(fname)
|
||||
|
||||
surf2 = ExplicitSurfaceThing(
|
||||
freq = 10,
|
||||
phase = 0,
|
||||
scale = 1/16, # from libfive
|
||||
inner = 0.9 * 1/16,
|
||||
outer = 2.0 * 1/16,
|
||||
rad = 0.3 * 1/16,
|
||||
ext_phase = numpy.pi/2)
|
||||
|
||||
z_init = 0
|
||||
x_init = surf2.x_cross(z_init, 1)
|
||||
mesh2 = surf2.screw_360(-5, x_init, z_init, 0.000002, 0.001, 5e-4)
|
||||
fname = "spiral_outer90_one_period.stl"
|
||||
print("Writing {}...".format(fname))
|
||||
mesh2.save(fname)
|
||||
Loading…
x
Reference in New Issue
Block a user