Welcome to Fractal Forums

Fractal Math, Chaos Theory & Research => Theory => Topic started by: Furan on February 22, 2013, 01:22:45 AM




Title: A generalization of triplex z^p+c
Post by: Furan on February 22, 2013, 01:22:45 AM
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

      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


Title: Re: A generalization of triplex z^p+c
Post by: kram1032 on February 22, 2013, 05:06:59 PM
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)


Title: Re: A generalization of triplex z^p+c
Post by: DarkBeam on February 22, 2013, 05:13:25 PM
this language is bizarre ;D heh


Title: Re: A generalization of triplex z^p+c
Post by: Furan on February 22, 2013, 09:55:15 PM
http://www.fractalforums.com/images-showcase-%28rate-my-fractal%29/furan%27s-corner-t13768/


Title: Re: A generalization of triplex z^p+c
Post by: kram1032 on February 23, 2013, 02:00:35 AM
(http://furan.sweb.cz/Wall2013/TrifoldMandelKnot960.jpg)?

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.


Title: Re: A generalization of triplex z^p+c
Post by: Furan on February 23, 2013, 01:42:24 PM
Good call, this does look a bit more interesting:
(http://furan.sweb.cz/Wall2013/TrifoldJuliaKnot.jpg)
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?


Title: Re: A generalization of triplex z^p+c
Post by: kram1032 on February 23, 2013, 05:57:15 PM
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)


Title: Re: A generalization of triplex z^p+c
Post by: Furan on February 23, 2013, 07:42:22 PM
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:

(http://furan.sweb.cz/Wall2011/GyroidCore01_400.jpg) (http://furan.sweb.cz/Wall2011/BlueGlassGyroid400.jpg)

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.


Title: Re: A generalization of triplex z^p+c
Post by: Kabuto on March 03, 2013, 12:13:22 AM
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...


Title: Re: A generalization of triplex z^p+c
Post by: kram1032 on March 03, 2013, 11:51:37 AM
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....)


Title: Re: A generalization of triplex z^p+c
Post by: Furan on March 03, 2013, 10:09:02 PM
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.


Title: Re: A generalization of triplex z^p+c
Post by: fractower on March 03, 2013, 10:36:56 PM
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?


Title: Re: A generalization of triplex z^p+c
Post by: kram1032 on March 04, 2013, 01:01:17 AM
objects of genus 2 or higher typically are split up into topological underwear (SFW) (http://en.wikipedia.org/wiki/Pair_of_pants). Using forms of that, you can generate any form you can imagine (and even those you can't), at least topologically.


Title: Re: A generalization of triplex z^p+c
Post by: DarkBeam on March 04, 2013, 09:40:13 AM
Can you rewrite without any use of GOTO, that's a bad practice of programming :) ... And in a single clear function? :)


Title: Re: A generalization of triplex z^p+c
Post by: Kabuto on March 04, 2013, 09:42:55 PM
@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?


Title: Re: A generalization of triplex z^p+c
Post by: kram1032 on March 05, 2013, 12:30:21 AM
Usually, Atan2(y,x) gives the angle as you want it. I find it highly unintuitive that it's not Atan2(x,y), but that's how, afaik, it's implemented in most languages.


Title: Re: A generalization of triplex z^p+c
Post by: Furan on March 05, 2013, 08:53:10 AM
There are actually several other errors in the code that didn't show up. I have to finish my presentation for tomorrow, than I'll be free to finish this. I corrected some errors but the outer cycle still diverges during first iteration. I'll have to check it with my initial algorithm in matlab.


Title: Re: A generalization of triplex z^p+c
Post by: Furan on March 08, 2013, 12:23:46 PM
Ok, its corrected now. Still not much difference. I was hoping a periodic pattern would emerge as in the MandelDonut.
Thank you, Kabuto, for pointing the error.
I'll try to improve the efficiency of the algorithm and put it here for others.
(http://furan.sweb.cz/Wall2013/TrifoldJuliaKnot1080.jpg)


Title: Re: A generalization of triplex z^p+c
Post by: kram1032 on March 08, 2013, 04:38:32 PM
But it IS periodic. Just not organized in bulbs as much as roots or vessels or something.
I like it very much! - and I'm curious what kind of crazy cross-combinations people will come up with once you got some correct and readible code :)


Title: Re: A generalization of triplex z^p+c
Post by: Alef on March 09, 2013, 01:41:44 PM

* 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)?
I think atan2 should computa angle between complex parts. That is complex argument.


Title: Re: A generalization of triplex z^p+c
Post by: kram1032 on March 09, 2013, 07:12:14 PM
yes, I believe, the question was about argument order.
Typically, if you have your complex plain and want your angle, beginning from +1, going to +i, then +1 and -i, Atan2 takes first the imaginary part and then the real one, so Atan2(y,x).