# Isosurface Dabbling

Following:
http://maxima.sourceforge.net/docs/manual/maxima_18.html
http://gkerns.people.ysu.edu/maxima/

load("vect");
load("eigen");

## Summary

I wanted to find a distance bound for my beloved 2005 isosurface that I used
in POV-Ray to make the corkscrew thingy. (Why: because I wanted to render
it directly in a sphere tracer in GLSL, or in Appleseed.)

I attempted to follow the method described in Hart's Sphere Tracing
paper, in which any implicit surface whose function is Lipschitz continuous
can be turned to a distance bound function for that same surface. This proved
very difficult since I am fairly certain that the function I created is *not*
Lipschitz continuous. There is also a pile of other failed methods here
after I dove down a rabbit-hole of vector calculus that I didn't really
understand very well, and then attempted to use CGAL to convert
this isosurface to a mesh, but it was creating meshes that were too
full of facetization artifacts (especially along the 'edge' of the corkscrew),
so I was trying to refine things. This is an unfinished mess, but the work
is left here.

At some point I realized my function already *is* an SDF. See section 2.1.
The formula works - I tried it in a sphere tracer in ShaderToy.
See: https://www.shadertoy.com/view/ttsXD4

The sections after 2.1 try to match up the 2005 formula and hand-derived
formula and prove that they are not just the same surface, but identical
formulas - provided that parameters are properly matched up.

Section 3.5 was my attempt to turn the implicit curve into an explicit
parametric one by using the cross-section. Section 3.6 was my attempt
to greatly simplify a bunch of annoying boundary conditions by using
an implicit curve for the XZ cross-section, and then just sweeping this
cross-section through the appropriate area.

### Questions

- How do I have an SDF that is not Lipschitz continuous?
Is that a 'sufficient, but not necessary' condition to produce a distance
estimate for a surface?

### Probably just me whining

- Can I do this in org-mode & babel? Jupyter? wxMaxima annoys me.
- Can I easily write this to let me work in x,y,z rather than
subscripts of a vector argument?

## Spiral Isosurface Thingy

as first explored in my 2005 notes

f3(p) := sqrt((p[2] * O - I * sin(p[1] * F + P))^2 + (p[3] * O - I * cos(p[1] * F + P))^2);

### SDFs derived by hand

I started from:
http://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm
and looked at opTwist & sdCylinder, which I was able to combine
to produce exactly the same surface.
The below is an exact SDF for cylinder along z (not a distance estimate).

cylSdf(p) := sqrt(p[1]^2+p[2]^2) - R;

and here is a domain transformation for a twist, starting with an XY rotation matrix
I use F2 instead of F.

rot : ev(matrix([C,-S, 0],[S,C,0],[0,0,1]), C=cos(F2*a-P), S=sin(F2*a-P));

twist(p) := ev(rot, a=p[3]) . [p[1], p[2], p[3]];

twist([x,y,z]);

so stacking a few transformations should produce an identical surface...

f : cylSdf(twist([x,y,z])  - [x0,y0,0])[1];

trigsimp(f);

This, tested in ShaderToy, works perfectly.
The constants may need some correction to match exactly what is used in my
2005 formula, but that should not be much of a problem.

### Why are these different?

The below is my 2005 formula. Only modification is that
I pass z,y,x instead of x,y,z.

g : ev(f3([z,y,x])-T);

trigsimp(g);

Evaluated at 4 phases of its full period:

makelist(ev(((g = 0) + T)^2,[z=i*%pi/(2*F)]), i, 0, 3);

and the derived SDF:

makelist(ev(((f = 0) + R)^2,[z=i*%pi/(2*F2)]), i, 0, 3);

Sorry for atrocious variable names - I was 18. See the following replacements:

replace : [y0=0, x0=I/O, R=(T/O), F2=-F];

...applied to the derived SDF:

trigsimp(ev(((f = 0) + R)^2, replace))*O^2;

This is indeed identical to the 2005 one:

sdf : ((trigsimp(g)=0)+T)^2;

### Analytical gradient

(still needs normalized. ugly denominators can probably be ignored.)

ev(express(grad(lhs(sdf))), diff);

gr : trigsimp(ev(express(grad(f)), diff));

optimize(gr);

## I suck at vector calculus

This was moved later on in the file on account of being
basically useless.

The point of this work is mostly to verify that some scratch
work I did with the CGAL library based around the code at
https://doc.cgal.org/latest/Surface_mesher/index.html#title3
is valid.
The issue I was trying to tackle is that the meshes it generated
were not quite optimal for me (perhaps due to me having different
goals than a package for computational geometry for CAD). The
positions of vertices in the mesh seemed to have a lot of random
noise, even if the overall error was bounded and could
be reduced arbitrarily by using more facets.  I was looking at
optimizing the mesh by doing per-vertex gradient descent to
nudge each vertex closer to the surface itself, which the
definition of an isosurface makes straightforward: for point
F(x,y,z), the closer F is to 0, the closer it is to the surface.
Gradient descent of course requires the gradient. I was dealing
with arbitrary chains of R^3 -> R^3 functions being composed with
an R^3 -> R function to produce the isosurface R^3 -> R, and didn't
feel like working out derivatives by hand.
Chris Hodapp, 2018-07-27

Note that the implicit surface itself is f3(x,y,z)=T where T is
some threshold. I use the squared version below because it doesn't
actually change the isosurface's properties, but is easier to work with.

surface : f3([x,y,z])^2 - T^2;

f_x : diff(surface, x);

f_y : diff(surface, y);

f_z : diff(surface, z);

f_yy : diff(f_y, y);

f_zz : diff(f_z, z);

f_yz : diff(f_y, z);

because I might forget:

f_zy : f_yz;

note that K below is the 2D curvature within the slice.,
it is not the 3D curvature of the surface!

K : (-f_y^2 * f_zz + 2*f_z*f_y*f_zy - f_z^2*f_yy) / ((f_z^2 + f_y^2)^(3/2));

K_y : diff(K, y);

K_z : diff(K, z);

we can try solving for singular points of the entire surface,
but while we can find singular points on some surface, and
can find points on the surface we want, we can't seem to find
both:

solve([f_x = 0, f_y = 0, f_z = 0], [x,y,z]);

solve([surface = 0], [x,y,z]);

solve([f_x=0, f_y=0, f_z=0, surface=0], [x,y,z]);

no singular points? may need lagrange multipliers after all...

guess I should start with the simpler 2D case.
that is: finding extrema of K under the constraint that surface=0.

extrema : solve([surface = 0, K_y = -L * f_y, K_z = -L * f_z], [L, y, z]);

length(extrema);

since we don't need the value of L:

[extrema[1][2], extrema[1][3]];

[extrema[2][2], extrema[2][3]];

[extrema[3][2], extrema[3][3]];

[extrema[4][2], extrema[4][3]];

those don't look bad.

params1 : [O=2.0, I=0.4, F=20, P=%pi, T=0.3];

surf1 : ev(surface, params1);

grad1 : ev([f_x, f_y, f_z], params1);

K1 : ev(K, params1);

wxplot3d(ev(K1, x=11), [y, -0.5, 0.5], [z, -0.5, 0.5],
        [zlabel,"K"]);

load(implicit_plot)$
wximplicit_plot(ev(surf1, x=0.22), [y, -0.5, 0.5], [z, -0.5, 0.5])$

pt : [x=-1.805000, y=-0.049920, z=0.005662];

ev(surf1, pt);

so, looks very close to surface

ev(grad1, pt);

so, nowhere near a singular point.

ev(K1, pt);

test1 : ev([extrema[3][2], extrema[3][3]], params1, pt[1]);

test2 : ev([extrema[4][2], extrema[4][3]], params1, pt[1]);

### wtf seriously

so, take a look again at surf1...

surf1 = 0;

IT'S A FREAKING CIRCLE

circle_params : [y2 = y - I*sin(P+x*F)/ O, z2 = z - I*cos(P+x*F) / O, R=T/O];

ev(y2^2 + z2^2 - R^2 = 0, circle_params, params1, eval);

and note that curvature is just 1/R:

ev(1/R, circle_params, params1, eval);

so the cross-section (sliced through x) is always a circle.
the 'sharpness' seen here has to be from something else, then.

what about if I use the full 3D curvature, simply using x as the tangent vector?

### 3D curvature

H_F : hessian(surface, [x, y, z]);

gradF : ev(express(grad(surface)), diff);

v: [vx, vy, vz];

Finally, normal curvature in 3D (for unit tangent vector v):

Kn : transpose(v) . H_F . v / sqrt(gradF . gradF);

Note that this is *not* the same as the curvature of any cross-section.
It is the curvature of: the curve projected onto the plane containing
the curve's tangent, and the surface normal.
(You don't choose the normal!)

Somewhat arbitrarily, I pick as my tangent vector: the cross product
of the surface normal (which is just grad F), and the direction the
cross-section moves (which is itself a tangent vector).

but how do I handle the v=[vx,vy,vz] here?
I need v declared as a vector but still something
I can evaluate like below:

ev(Kn, vx=v_xc[1], vy=v_xc[2], vz=v_xc[3], params1, eval);

--- ignore the below, it's earlier scratch work ---

is this the tangent I want? I'm pretty sure it's not since, you know,
IT'S NOT AN ACTUAL TANGENT OF THE SURFACE.
needs to be perpendicular to the normal, doofus.

Kx : f_xx / (f_x^2 + f_y^2 + f_z^2)^0.5;

Kx_x : diff(Kx, x)$
Kx_y : diff(Kx, y)$
Kx_z : diff(Kx, z)$

the below hangs for some reason (possibly because my tangent vector is completely wrong):

extrema2 : solve([surface = 0, Kx_y = -L * f_y, Kx_z = -L * f_z], [L, y, z]);

### Lipschitz constant or something

So, the gradient is unbounded in each component.
Is this the same as the total derivative being unbounded?
If that's the case, then the function isn't Lipschitz continuous.
But... various other functions express the same surface. Like this:

surf2 : f3([x,y,z]) - T;

ev(express(grad(surf2)), diff);

This has order 1 on top, order 1 on bottom, so this should balance out.
If I understand right what Wikipedia is saying with,
"Moreover, if K is the best Lipschitz constant of f, then |Df(x)| ≤ K
whenever the total derivative Df exists."
then...

grad2 : ev(express(grad(surf2)), diff)$

Df : trigsimp(sqrt(grad2 . grad2));

ev(Df, x=-P/F);

trigsimp(ev(Df, x=-P/F + %pi / F / 2));

trigsimp(ev(Df, x=-P/F + %pi / F));

trigsimp(ev(Df, x=-P/F + 3 * %pi / F / 2));

this is guesswork under the assumption that max/min occur where sin/cos
hit their max/min

Df2 : subst([sin(P+x*F)=S, cos(P+x*F)=C], Df);

So I *think* that regardless the additive terms in the numerator -
(z^2-y^2)*F^2*I^2*O^2, 2*y*z*F^2*I^2*O^2, (z^2+y^2)*O^4, (y^2*F^2+1)*I^2*O^2, and  -
will dominate due to the powers on y and z.

Denominator is dominated by (-z^2-y^2)*O^2.

(z^2-y^2)*F^2*I^2*O^2 becomes (for large z or y) I^2 * F^2.
2*y*z*F^2*I^2*O^2 vanishes? I think?
(z^2+y^2)*O^4 becomes (for large z or y) O^2.
(y^2*F^2+1)*I^2*O^2 becomes (for large y) I^2 * F^2.

Then, if I'm right, the overall bound is something like: |O|+|I*F|
This seems sensible, I think.
Intuitively: higher (spatial) frequencies correspond to steeper
change.

ev(Df2, [z=0, y=0]);

ev(Df2, [z=1, y=0]);

ev(Df2, [z=0, y=1]);

### More Lipschitz derivation?

The below is a lower-bound on K:

K : ((f3([x1,y1,z1]) - T) - (f3([x2,y2,z2]) - T)) / sqrt([x1,y1,z1] . [x2,y2,z2]);

trigsimp(K);

So is this (it's just squaring numerator & denominator) but I feel like
I just made it worse:

K2 : ((f3([x1,y1,z1]) - T) - (f3([x2,y2,z2]) - T))^2 / ([x1,y1,z1] . [x2,y2,z2]);

trigsimp(K2);

If one considers that cos/sin are both bounded to [-1,1]...

K_simp : subst([sin(P+x1*F)=S1, cos(P+x1*F)=C1, sin(P+x2*F)=S2, cos(P+x2*F)=C2], trigsimp(K));

Also, the trig terms are always paired with a lower-order y or z
(they aren't squared, but are under a radical) than those of the
denominator.

trigsimp(ev(K, [y1=0, y2=0, z1=0, z2=0]));

so, this is at a minimum for those...
but suppose we let (for instance) z1 or z2 grow without bound.

trigsimp(ev(K, [y1=0, y2=0, z1=0, z2=inf]));

trigsimp(ev(K, [y1=0, y2=0, z1=inf, z2=0]));

then K grows without bound (nothing in the numerator counteracts it,
even if the other z is nonzero). 

trigsimp(ev(K, [y1=0, y2=0, z1=1, z2=inf]));

but... as long as f^-1(x) remains the same, this is the same surface.
I can just take it to a lower power.

f4(p) := ((p[2] * O - I * sin(p[1] * F + P))^2 + (p[3] * O - I * cos(p[1] * F + P))^2)^(1/4);

K4 : ((f4([x1,y1,z1]) - T) - (f4([x2,y2,z2]) - T)) / sqrt([x1,y1,z1] . [x2,y2,z2]);

trigsimp(ev(K4, [y1=0, y2=0, z1=inf, z2=0]));

oh, I guess that doesn't really do much since the infinity is still there
in the numerator.....

so it is *not* Lipschitz continuous, even if its isosurface is continuous.
It is probably locally Lipschitz continuous but I don't know that that
gets me much further than simply ray-marching.

### Other cross sections

This is for y = 0.
Imagine spinning the entire graph around the x axis, and as you
do, moving it to the right (positive x) so that a half-turn
makes each 'bottom' lobe line up with the next 'top' lobe to the
right, and a full turn makes the 'bottom' lobe line up with the
next bottom one (and the top one with the next top one).

wximplicit_plot(ev(surf1, y=0), [x, -0.5, 0.5], [z, -0.5, 0.5])$

I could probably have inverted this in Maxima, but instead did it by hand:

x_inv : acos((T^2-I^2-(O*z)^2)/(2*O*z*I))/F - P/F;

arccos is periodic with integer multiples of 2*pi, so I added that below:

wxplot2d([ev(x_inv, params1), ev(x_inv + 2*%pi/F, params1), ev(x_inv - 2*%pi/F, params1)], [z, -0.5, 0.5])$

This looks right except that I'm losing an entire half of each lobe and I can't
reason through why.
If cos(x)=y, then x=arccos(y)+2*pi*n.
But: if cos(x)=y then cos(-x)=y, so y=-arccos(x)+2*pi*n.
So, below I have negated F and P (which is equivalent to negating (F*x+P),
the argument to arccos:

x_inv2 : ev(x_inv, F=-F, P=-P);

wxplot2d([ev(x_inv2, params1), ev(x_inv2 + 2*%pi/F, params1), ev(x_inv2 - 2*%pi/F, params1)], [z, -0.5, 0.5])$

so there we go:
(how can I make wxplot2d get rid of the annoying labels?)

wxplot2d([ev(x_inv, params1), ev(x_inv2, params1)], [z, -0.5, 0.5])$

z_inv : solve(ev(surface, y=0), z);

oh. duh - it was a quadratic, and I ignored it the first time because
I saw that it wouldn't factor easily.
well, regardless, it has to be piecewise even for just one period worth,
like the other.

z_inv_ev : [rhs(ev(z_inv, params1)[1]), rhs(ev(z_inv, params1)[2])];

wxplot2d(z_inv_ev, [x, -0.5, 0.5])$

For the sake of my explicit representation to generate a mesh, I need to guarantee that a
vertex is directly at the valley at the bottom, and somewhat also at the peak on top.

solve(rhs(diff(z_inv[1], x)) = 0, x);

solve(rhs(diff(z_inv[2], x)) = 0, x);

this first solution checks out: if phase P=0, x=0 is a peak of both.

I may have one additional issue: the vertex sampling needs to be sensible,
as if I should be incrementing x *inversely* proportionally with the
derivative of f(x).
At first glance, my vertices should be placed at (approximately) equal
arc lengths along the function itself, but this also isn't fully right,
because I also want portions with higher curvature to have more vertices
devoted to them.

Derivative:

wxplot2d(diff(z_inv_ev, x), [x, -0.5, 0.5])$

Curvature:

Kz : diff(z_inv, x, 2) / (1 + diff(z_inv, x)^2)^(3/2);

optimize(Kz);

Kz : diff(z_inv_ev, x, 2) / (1 + diff(z_inv_ev, x)^2)^(3/2)$
wxplot2d(Kz, [x, -0.05, 0.2])$

ev(Kz, x=0);

Arc length:

Sz1(b) := quad_qag(sqrt(1+diff(z_inv_ev[1],x)^2), x, 0, b, 3)[1];

Sz2(b) := quad_qag(sqrt(1+diff(z_inv_ev[2],x)^2), x, 0, b, 3)[1];

wxplot2d([Sz1(b), Sz2(b)], [b, 0, 0.04])$

I guess that's nothing particularly interesting.

### Ugh, just use the implicit curve

I switched x and z to match my code. Sorry, not sorry.

xz_imp : ev(surface, [y=0, x=z, z=x]);

F_x : trigsimp(diff(xz_imp, x));

F_z : trigsimp(diff(xz_imp, z));

F_xx : trigsimp(diff(F_x, x));

F_zz : trigsimp(diff(F_z, z));

F_xz : trigsimp(diff(F_x, z));

K : trigsimp((-F_z^2*F_xx+2*F_x*F_z*F_xz-F_x^2*F_zz)/((F_x^2+F_z^2)^(3/2)));

for ex in optimize(K) do print(lhs(ex),"=",rhs(ex))$

Tentative approach:
- I have the gradient with F_x & F_z. Moving perpendicular to
the gradient should mean moving along the implicit curve. If I
keep the step size small, I should then be able to follow the
gradient again to adjust my point back to the implicit curve.
- I now have a non-piecewise K and can adjust my step size
accordingly.
- I need a seed point (where F(x,z)=0) but that's pretty easy.

## where I'm stuck and/or what I learned

- Cross section curvature was completely the wrong way to go here
- 3D curvature has its own complexities
- I can't see how to solve for that d vector (direction of the cross section)
without finding dy/dx and dz/dx of F, which is tricky when it's an implicit
function, though this may still be possible
- I already know that while the edges are sharp, they are neither 1D features
nor singular points/curves and it may not even have the desired effect to
make CGAL preserve them
- I am pretty sure I can find an analytical answer without requiring
curvature

- Can't I find a way to make CGAL subdivide more in areas of higher curvature,
or something like that?  Could I subdivide it myself and then do
gradient-descent?