A generalization of triplex z^p+c

[Furan]

Let us imagine a parametric function X=f(r,phi,theta)
(A surface at a threshold r and two angles phi, theta.)

Let there be G = [r,theta,phi] = g(X) = f^-1(X), either explicit or implicit. In case f^-1 is not a function, g could be altered to return the solution with minimal r.

With possible adjustment of the constant C=c(X_0), iterate:
G     <-  f^-1(X)
X     <-  f(r^p,p*phi,p*theta) + C

Mandelbulb for f(G) as a sphere, MandelDonut for a torus. What else?

Trifold MandelKnot

The true Trifold Knot can be expressed probably only as a parametric function.
To get X from r,phi,theta, an initial guess is acquired via an empirical function.
Following Newton-Rhapson:
G_n+1 = G_n + J^-1 * (X-f(G_n)); where J is the Jacobian: partial df_i/dG_j

Previous step should work (sometimes) as initial guess. However something isn't working, I'll keep trying.
To get the idea:

Initial guess at parameters r,phi,theta

Code:

  100 phi = atan2(x,y)
C
      t2 =  cos(1.5*phi)/3.8
      t1 =  pi/12.*cos(1.5*phi-t2) - phi/2. + pi/2.
      t2 = -pi/12.*cos(1.5*phi+t2) - phi/2. + 1.5*pi
C
      X1=3.*(/sin(t1)+2.*sin(2.*t1),cos(t1)-2.*cos(2.*t1),-sin(3.*t1)/)
      X2=3.*(/sin(t2)+2.*sin(2.*t2),cos(t2)-2.*cos(2.*t2),-sin(3.*t2)/)
      qnX1=(x-X1(1))**2+(y-X1(2))**2+(z-X1(3))**2
      qnX2=(x-X2(1))**2+(y-X2(2))**2+(z-X2(3))**2

      IF (qnX1.LT.qnX2) THEN
        r = sqrt(qnX1)
        phi = t1
        theta=-atan2(z-X1(3),sqrt(x**2+y**2)-sqrt(X1(1)**2+X1(2)**2))
      ELSE
        r = sqrt(qnX2)
        phi = t2
        theta=-atan2(z-X2(3),sqrt(x**2+y**2)-sqrt(X2(1)**2+X2(2)**2))
      END IF
      j=1

Iterate to get parameters r,phi,theta

Code:

  110 sphi(1)=sin(phi)          !sin/cos lookup table
      sphi(2)=sin(2.*phi)
      sphi(3)=sin(3.*phi)
cphi(1)=cos(phi)
cphi(2)=cos(2.*phi)
cphi(3)=cos(3.*phi)
C
      v0=1./sqrt(32.*cphi(3)+68.) !inverted norms
      w0=1./sqrt(8.)/
     /   sqrt(795.+36.*cos(9.*phi)+217.*cos(6.*phi)+652.*cphi(3))
C
      v1=(-2.*sphi(1)+8.*sphi(2))*v0
      v2=-(2.*cphi(1)+8.*cphi(2))*v0
C
      w1=-12.*cphi(3)*(cphi(1)+4.*cphi(2))*w0
      w2= 12.*cphi(3)*(sphi(1)-4.*sphi(2))*w0
      w3= (-68.-32.*cphi(3))*w0
C
      cth=cos(theta)
      sth=sin(theta)
C     New iteration
      Xn   = (/a*sphi(1)+2*a*sphi(2)+r*(v1*cth+w1*sth),
     ,         a*cphi(1)-2*a*cphi(2)+r*(v2*cth+w2*sth),
     ,        -a*sphi(3)+r*w3*sth/)
dX   = (/Xn(1)-x,Xn(2)-y,Xn(3)-z/)
      qndX = dX(1)**2 + dX(2)**2 + dX(3)**2
j=j+1
IF ((qndX.LT.1.D-24).OR.(j.LT.30)) GOTO 113

Calculating Jacobian + some tweaks to get it stable.

Code:

C      dX/dr
      Jac(1,1) = v1*cth+w1*sth
      Jac(2,1) = v2*cth+w2*sth
      Jac(3,1) = w3*sth
C      dX/dtheta
      Jac(1,3) = -r*(v1*sth+w1*cth)
      Jac(2,3) = -r*(v2*sth+w2*cth)
      Jac(3,3) =  r*w3*cth
C
      dv0dphi = -v0*48*sphi(3)
      dv1dphi = (-2*cphi(1)+16*cphi(2))*v0 - v1*v0*dv0dphi
      dv2dphi = ( 2*sphi(1)+16*sphi(2))*v0 - v2*v0*dv0dphi
C 
      dw0dphi = 4*w0*(-324*sin(9*phi)-1302*sin(6*phi)-1956*sin(3*phi))
      dw1dphi = (36*sphi(3)*(cphi(1)+4*cphi(2))-12*cphi(3)*
     *          (-sphi(1)-8*sphi(2)))*w0-w1*dw0dphi*w0
      dw2dphi = (-36*sphi(3)*(sphi(1)-4*sphi(2))+12*cphi(3)*
     *          (cphi(1)-8*cphi(2)))*w0-w2*dw0dphi*w0
      dw3dphi = 96*sphi(3)*w0-w3*w0*dw0dphi
C      dX/dphi
      Jac(1,2) =  a*cphi(1)+4*a*cphi(2)+r*cth*dv1dphi+r*sth*dw1dphi
      Jac(2,2) = -a*sphi(1)+4*a*sphi(2)+r*cth*dv2dphi+r*sth*dw2dphi
      Jac(3,2) = -3*a*cphi(3)+r*sth*dw3dphi
C
      invJac(1,1) = Jac(2,2)*Jac(3,3)-Jac(3,2)*Jac(2,3)
      invJac(1,2) = Jac(3,2)*Jac(1,3)-Jac(1,2)*Jac(3,3)
      invJac(1,3) = Jac(1,2)*Jac(2,3)-Jac(2,2)*Jac(1,3)
      detJac = 1./(Jac(1,1)*invJac(1,1)+Jac(2,1)*invJac(1,2)+
     +             Jac(3,1)*invJac(1,3))
C Inverse Jacobian
      invJac(1,1) = invJac(1,1)*detJac
      invJac(1,2) = invJac(1,2)*detJac
      invJac(1,3) = invJac(1,3)*detJac
      invJac(2,1) = (Jac(2,3)*Jac(3,1)-Jac(2,1)*Jac(3,3))*detJac
      invJac(2,2) = (Jac(1,1)*Jac(3,3)-Jac(1,3)*Jac(3,1))*detJac
      invJac(2,3) = (Jac(1,3)*Jac(2,1)-Jac(1,1)*Jac(2,3))*detJac
invJac(3,1) = (Jac(2,1)*Jac(3,2)-Jac(2,2)*Jac(3,1))*detJac
      invJac(3,2) = (Jac(1,2)*Jac(3,1)-Jac(1,1)*Jac(3,2))*detJac
      invJac(3,3) = (Jac(1,1)*Jac(2,2)-Jac(1,2)*Jac(2,1))*detJac
      dr     = invJac(1,1)*dX(1)+invJac(1,2)*dX(2)+
     +         invJac(1,3)*dX(3)
dphi   = invJac(2,1)*dX(1)+invJac(2,2)*dX(2)+
     +         invJac(2,3)*dX(3)
dtheta = invJac(3,1)*dX(1)+invJac(3,2)*dX(2)+
     +         invJac(3,3)*dX(3)
IF (dtheta.LT.0.1) dr=0.5*dr*cos(dtheta)
      IF (abs(dtheta).GT.1.5) dtheta=1.5*sign(1.,dtheta)
      theta=mod(theta+dtheta,2*pi)
      IF (dr+r.LT.0.) dr=0
      r=r+dr
      IF (r.GT.25) r=25
      IF (abs(dphi).GT.0.1) dphi=0.1*sign(1.,dphi)
      phi=phi+dphi
      GOTO 110

And the very fractal

Code:

  113 IF (ni.EQ.1) THEN
        c=r
END IF
        xc=r*(v1*cth+w1*sth)
        yc=r*(v2*cth+w2*sth)
        zc=r*w3*sth
      r=r**3
IF ((r.GT.10.).OR.(ni.EQ.nimax)) GOTO 115
      phi=phi*3.
theta=theta*3.
sphi(1)=sin(phi)          !sin/cos lookup table
      sphi(2)=sin(2.*phi)
      sphi(3)=sin(3.*phi)
cphi(1)=cos(phi)
cphi(2)=cos(2.*phi)
cphi(3)=cos(3.*phi)
C
      v0=1./sqrt(32.*cphi(3)+68.) !inverted norms
      w0=1./sqrt(8.)/
     /   sqrt(795.+36.*cos(9.*phi)+217.*cos(6.*phi)+652.*cphi(3))
C
      v1=(-2.*sphi(1)+8.*sphi(2))*v0
      v2=-(2.*cphi(1)+8.*cphi(2))*v0
C
      w1=-12.*cphi(3)*(cphi(1)+4.*cphi(2))*w0
      w2= 12.*cphi(3)*(sphi(1)-4.*sphi(2))*w0
      w3= (-68.-32.*cphi(3))*w0
C
      cth=cos(theta)
      sth=sin(theta)
C     New fractal iteration
      x = a*sphi(1)+2*a*sphi(2)+r*(v1*cth+w1*sth) + xc
      y = a*cphi(1)-2*a*cphi(2)+r*(v2*cth+w2*sth) + yc
      z = -a*sphi(3)+r*w3*sth                     + zc
C
ni=ni+1
GOTO 100

[kram1032]

You forgot the most important part. How does it look like?

(Notably, if you had given images but no implementation details, you would also have forgotten the most important part)

[DarkBeam]

this language is bizarre heh

[Furan]

http://www.fractalforums.com/images-showcase-%28rate-my-fractal%29/furan%27s-corner-t13768/

[kram1032]

?

Looks really nice!
Do you think you could generalize the very same thing to tensor product surfaces of any shape? - in that case you could use a huge number of generic shapes as base for Mandelbulb-like attractors!

Also, I'd love to see Juliasets of this.

[Furan]

Good call, this does look a bit more interesting:

Julia: C = [0.5 0.5 0.5] The bigger the C, the more sparse is the object.
I will run a high-res version in a few days. I'm still working on the coloring scheme and optimization.

This is the first time I heard about tensor product surfaces. Do you have a favorite one I could try to implement?

[kram1032]

My understanding of this is a bit limited so I might confuse stuff here and it won't work as I imagine, but:
Tensor-product-surfaces can be used as a form of free-form surfaces. If I'm not mistaken, you could essentially use them to get any polygonal mesh into a Mandelbulb-esque form.

(I really hope I'm not mixing apples and oranges here)

[Furan]

Well, even when divided into spline patches, the surface would have to have a parametric representation (analytical, or numerical).
That would be probably impossible for objects like these two:



Or another approach: For a point in space X, the original implicit object f(X)=0 could return grad f and "eject" the point outside in the opposite direction according to the value of f(X). It could work as long as the ejected distance function formed an attractor. This would not be fractal-like in every point.

[Kabuto]

Furan: Nice fractal, but I don't understand your code. According to it it consists of 2 loops: an outer one (fractal iterations) and an inner one (to optimize parameters). However, the following line of code:

IF ((qndX.LT.1.D-24).OR.(j.LT.30)) GOTO 113

implies that the inner loop is left on its first iteration and that thus the entire Jacobian part is dead code...

[kram1032]

what the -- I did post a reply here days ago and only now I see that it's gone.

Those two glass objects look awesome too. Couldn't you approximate parametric coordinates for stuff like that?
Or choose "arbitrary" parametric coordinates like you'd do for textures? (E.g. you could use UV coordinates by unwrapping the mesh? - possibly as auto-generated by something like ptex? I'm really not sure how you'd go about finding inverses and such for that then though....)

[Furan]

Kabuto: I can't believe this. I remember checking inner iteration count. I probably made a mistake and didn't even realize it. I'll check it. I'll be back tomorrow.

kram1032: As I was saying, I don't think it's possible. And even if you could create a mesh and unwrap it, it would be discontinuous and very wildly distorted.

[fractower]

What are the topological limits of closed surface parameterization. I suppose the sphere and torus are easy. I have seen Klein bottle. Can a two hole closed object be done with two variables or is that the breaking point?

[kram1032]

objects of genus 2 or higher typically are split up into topological underwear (SFW). Using forms of that, you can generate any form you can imagine (and even those you can't), at least topologically.

[DarkBeam]

Can you rewrite without any use of GOTO, that's a bad practice of programming ... And in a single clear function?

[Kabuto]

@Furan: are there any specs for this pretty uncommon BASIC-like language? I've never seen this syntax before.

Some specific questions where I can only make guesses:

* What exactly does atan2 do? Does atan2(y,x) behave like atan(y/x)+(x < 0 ? pi*sign(y) : 0) for x!=0 and y!=0? Or is it rather the other way round, i.e. atan2(x,y)?
* Where does a come from? (Used after the comment line "New iteration")
* Does dr+r.LT.0. actually mean (dr+r) < 0 ?
* Is there any reason behind constants being floats (e.g. 4. or 10.) in some parts and integers (e.g. 4 or 10) in other parts?