Logo by Fiery - Contribute your own Logo!

END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

it was a great time but no longer maintainable by c.Kleinhuis contact him for any data retrieval,
thanks and see you perhaps in 10 years again

this forum will stay online for reference
News: Check out the originating "3d Mandelbulb" thread here
 
*
Welcome, Guest. Please login or register. November 19, 2018, 07:40:24 PM


Login with username, password and session length


The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!


Pages: [1] 2   Go Down
  Print  
Share this topic on DiggShare this topic on FacebookShare this topic on GoogleShare this topic on RedditShare this topic on StumbleUponShare this topic on Twitter
Author Topic: smooth external angle of Mandelbrot set?  (Read 4679 times)
0 Members and 1 Guest are viewing this topic.
billtavis
Safarist
******
Posts: 96


WWW
« on: May 06, 2014, 12:07:05 AM »

Does anyone know of a way to calculate a smooth field representing the external angle of the Mandelbrot set? Or any source code to point me towards? I have found many mathematical explanations but a lot of it is over my head.  huh?  I'm not looking for the fastest way, just the easiest to understand.

This site http://eventoneresearch.com/ has exactly what I'm looking for, the "field gradient". But it looks like their technique is patented. Any help is much appreciated, thanks!
Logged
youhn
Fractal Molossus
**
Posts: 696


Shapes only exists in our heads.


« Reply #1 on: May 06, 2014, 07:51:34 AM »

Check out Gnofract4D. Has lots of coloring methods, all open source.
Logged
hobold
Fractal Bachius
*
Posts: 573


« Reply #2 on: May 06, 2014, 09:54:52 AM »

I thought that should be simple, based on a "binary decomposition". Every iteration provides another bit of external angle, from most significant to least significant. It doesn't quite work, but it comes tantalizingly close. Here is the code anyway, in case I only made some stupid error which others might spot.

Code:
// Mandelbrot iteration of a single point
// tries to compute angle of external ray passing through that point (buggy!)
unsigned iterate(const double cx, const double cy, const long maxiter) {
    unsigned i = 0;                 // iteration count
    double zx = 0.0, zy = 0.0;      // complex z
    double z2x = 0.0, z2y = 0.0;    // z real and imaginary parts squared
   
    unsigned angle = 0;             // accumulates fractional angle
    unsigned curBit = 0x80000000UL; // current bit of binary fraction
   
    do {
        i++;
        if (i > maxiter) {
            return -1;      // alternative loop end
        }
       
        z2x = zx*zx;        // squared components
        z2y = zy*zy;

        zy = 2.0*zx*zy;     // complex square, in place
        zx = z2x - z2y;
       
        zx = zx + cx;       // plus constant
        zy = zy + cy;
       
        // accumulate another bit of fractional angle
        if (zy < 0.0) {
            angle |= curBit;
        }
        curBit = curBit >> 1;   // next bit
       
    } while ((zx*zx + zy*zy) < 10000.0);
   
    // keep refining the angle for escaped points
    while (curBit > 0) {
        z2x = zx*zx;        // squared components
        z2y = zy*zy;
       
        zy = 2.0*zy*zx;     // complex square
        zx = z2x - z2y;
       
        zx += cx;           // plus constant
        zy += cy;
       
        // accumulate another bit of fractional angle
        if (zy < 0.0) {
            angle |= curBit;
        }
        curBit = curBit >> 1;   // next bit
    }
   
    return angle;
}

But it could also be the case that my understanding of the idea only really applies to Julia sets, not to the Mandelbrot set.
Logged
claude
Fractal Bachius
*
Posts: 563



WWW
« Reply #3 on: May 06, 2014, 06:29:11 PM »

I'm not looking for the fastest way, just the easiest to understand.

Here's what I did:

For each exterior pixel, trace its external ray through the C plane outwards accumulating bits in reverse order as you pass each dwell band - once the ray endpoint exceeds the bailout radius then reverse the bitstring.  That gives you a path through the binary decomposition grid, you can use the final angle of the original iterate (once it passes the bailout radius) to get local coordinates within each binary decomposition cell for a bit more smoothing.

http://www.math.nagoya-u.ac.jp/~kawahira/programs/mandel-exray.pdf has an algorithm for tracing rays from outwards in, it needs a bit of fiddling to trace from inwards out - in particular when crossing dwell bands you need to check which of the two angle doubling preimages is closer to the current ray endpoint (which you choose determines one bit of the output angle bitstring).

Here's one image I made (still from a video), using the angle field and potential to map a hyperbolic tiling to the exterior:
http://mathr.co.uk/mandelbrot/2013-08-20_mandelbrot_exterior_hyperbolic_tiling.png
Logged
laser blaster
Iterator
*
Posts: 178


« Reply #4 on: May 06, 2014, 07:51:32 PM »

Wow, that's really cool! This opens up the door for some cool texture effects.
Logged
billtavis
Safarist
******
Posts: 96


WWW
« Reply #5 on: May 06, 2014, 09:31:39 PM »

thanks everyone for the responses! I've got a great place to start from with the explanations and code
Logged
David Makin
Global Moderator
Fractal Senior
******
Posts: 2286



Makin' Magic Fractals
WWW
« Reply #6 on: June 12, 2014, 03:02:20 PM »

I realise I'm a little late to this but I've been away from fractaling for a while.
I wrote this a while ago and I spent a long time fiddling to get my version correct to as great a depth as possible.

Here's my code from mmf3.ucl for Ultra Fractal:

*************
Code:
MMF3-FieldLines {
;
; This version of calculating field lines will work reasonably well with both
; Mandelbrot and Julia Sets for divergent fractals with a divergence that is
; close to being a positive integer >=2 but can also produce quite interesting
; results with fractals that do not fit that criteria - just not rendering the
; field lines correctly in other cases :-)
;
; Thanks to Chris Hayton for showing me how to fix the "off-by-one" problems.
;
global:
  float twopi = 2.0*#pi
  float dtwopi = 0.5/#pi
  float dp = 1.0/@power
  float m = @power*dtwopi
  float dp2 = 1.0/(2.0+@power) ; Changed from originally trying 1/power
init:
  complex zv[#maxiter]
  int i = 0
  float a = 0.0
  float b = 0.0
  float c = 0.0
  float s = 0.0
  if @julia
    zv[i] = #z
    i = i + 1
  endif
loop:
  zv[i] = #z
  i = i + 1
final:
  if (s = atan2(#z))<0
    s = s + twopi
  endif
  while i > @skipiter
    i = i - 1
    if @version==0
      b = (@power*atan2(zv[i]))%twopi
    else
      b = atan2(zv[i]^@power)
    endif
    if b < 0.0
      b = b + twopi
    endif
    c = c + b - s
    if c >= #pi
      s = s + twopi
      if @fixit || @accident
        c = c - twopi
      endif
    elseif c < -#pi
      s = s - twopi
      if @fixit || @accident
        c = c + twopi
      endif
    endif
    if @fixit
      c = dp2*c
    elseif !@accident
      c = 0.0
    endif
    if (a = atan2(zv[i])) < 0.0
      a = a + twopi
    endif
    s = dp*(s + twopi*floor(a*m))
  endwhile
  s = s*dtwopi
  if @fixskip
    while s<0.0
      s = s + 1.0
    endwhile
    while s>=1.0
      s = s - 1.0
    endwhile
  endif
  #index = s
default:
  title = "MMF3 Field Lines (use high bailouts)"
  heading
    caption = "Information"
;    text = "This colouring is a little touchy around certain values, if you \
;            get obvious lines or dots that shouldn't be there but the 'Degree \
;            of Divergence' is set properly for the main formula and you are \
;            using a high bailout of say 1e20 or higher then try \
;            adjusting the location very slightly. For example if vertical \
;            lines appear on an unrotated fractal then try changing the x \
;            (real) location very slightly, if horizontal lines appear then try \
;            changing the y (imag) location very slightly."
  endheading
  heading
  endheading
  param version
    caption = "Version"
    enum = "Original" "Alternative"
    default = 1
    hint = "Allows you to modify one of the calculations in the formula that \
            is particularly sensitive to small errors. Ideally both methods \
            should produce identical results but they don't, so if you get \
            visible errors using one method then try the other (assuming you \
            have the divergence set correctly and are using a high bailout)."
  endparam
  param power
    caption = "Degree of divergence"
    default = 2.0
    hint = "This is an estimate of the divergence of the main formula, e.g. 2 \
            for z^2+c, 3 for z^3+c etc."
  endparam
  param julia
    caption = "Julia ?"
    default = false
    hint = "Enable for correct rendering of Julia Sets."
  endparam
  param skipiter
    caption = "'Start' iteration"
    default = 0
    min = 0
    hint = "You can use this to attempt to get better detail at higher \
            iteration depths, note that it only works as intended with 'Color \
            Density' as a whole integer >=1."
  endparam
  param fixskip
    caption = "Fix 'Start' iteration"
    default = false
    hint = "Fixes a problem when the start iteration is non-zero and maybe in \
            other cases - specifically if you get solid areas of palette \
            colour zero when the Transfer Function is Linear this may fix the \
            problem."
  endparam
  param fixit
    caption = "Fix Field Lines"
    default = false
    hint = "Enable this to make the formula much more accurate at rendering \
            the field lines correctly. It's a fudge produced after trial and \
            error but is quite effective at correcting some areas."
  endparam
  param accident
    caption = "Happy Accident"
    default = false
    hint = "You may find the colouring useful when this is enabled. This \
            is available as I found the results when initially adding the \
            'Fix it' option quite interesting."
    visible = !@fixit
  endparam
}

**************
Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
David Makin
Global Moderator
Fractal Senior
******
Posts: 2286



Makin' Magic Fractals
WWW
« Reply #7 on: June 12, 2014, 08:21:26 PM »

Theere's a later version of my formula for field lines - though somewhat extended to include smooth iteration, distance estimation and distance estimation angles in my class formulas for UF - mmf.ulb called "MMF Field Estimator"
It's a UF class formula so not quite so easy to follow but I think it should be clear enough.
I boiled down "corrections" to field lines that work in most places but at the moment require manual adjustment to work properly - these are the two parameters that can be modified by the user to fix errors in the field lines, unfortunately you can't always fix them everywhere in view even then but this is about as close as you can get I think.
I tried including the MMF Field Lines code here but unfortunately it's too long, it's in mmf.ulb in te Ultra Fractal Formula database at:

http://formulas.ultrafractal.com/
Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
billtavis
Safarist
******
Posts: 96


WWW
« Reply #8 on: June 13, 2014, 07:00:37 PM »

thanks David! I'll look through that code and see what I can get from it, what I have now isn't entirely satisfactory. I guess I should be learning more about UltraFractal as well! I don't think I've ever tried that one
Logged
billtavis
Safarist
******
Posts: 96


WWW
« Reply #9 on: February 19, 2015, 05:59:47 PM »

Thank you everyone for your input. I was finally able to get what I was looking for: a very smooth external angle rendering that I can understand (and yes it's very slow to calculate). I solved it by pretending that from every pixel, I drop a ball and let it roll downhill away from the set. Once it reaches a very large distance from the set, I use it's angle on that large circle. Every step of the ball rolling requires computing the gradient at its current location, which requires a full set of iterations. This is why it's so slow. However the advantage is that by decreasing the step size and using 4th order Runge-Kutta integration I was able to get it quite accurate:

Logged
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« Reply #10 on: February 19, 2015, 06:26:38 PM »

cool. Src code with extensive description is welcome
Logged
youhn
Fractal Molossus
**
Posts: 696


Shapes only exists in our heads.


« Reply #11 on: February 19, 2015, 06:27:32 PM »

That looks nice, like arms scanning/touching the 2D boundary of the set.  smiley

Also reminds me of synapsis, or how do you call those neuron-ends in our brain?
Logged
DarkBeam
Global Moderator
Fractal Senior
******
Posts: 2512


Fragments of the fractal -like the tip of it


« Reply #12 on: February 19, 2015, 06:47:23 PM »

Is this related to triangle inequality average? No... wink
Logged

No sweat, guardian of wisdom!
billtavis
Safarist
******
Posts: 96


WWW
« Reply #13 on: February 20, 2015, 05:48:27 PM »

I don't know about triangle inequality average, the pictures I found looks related but different.

So the first thing I did was to get a binary version working, based on adding bits to a fractional angle. I used the source code found here: http://www.mndynamics.com/indexp.html to get a working version. The resulting grayscale image was very smooth (shown by the contour lines) but it's pretty obvious that it doesn't capture the complexity of the boundary.


Then I found a formula for calculating the gradient of the potential (smoothed iteration count) at this site: http://linas.org/art-gallery/escape/ray.html  So, for a given c-value in the plane, I would get the distance to the set using the standard distance estimation, and then take a step of that length going down the gradient. This would give me a new c-value, where I would find the new distance and gradient and take another step, until c was very far from the origin. Then, I would plot the original pixel with the angle of the distant c-value. The resulting grayscale image looked just as smooth as the binary version, and it captured the detail at the boundary. However, when I ran contour lines along it there was obvious errors creeping in from taking too large of a step.


I realized that this was a case of Euler integration, so I knew I could improve the accuracy by switching to Runge-Kutta instead. Also, I limit the step size to a fraction of the distance, and I give a maximum cutoff. Limiting the step size along with using Runge-Kutta gave a result just as smooth as the binary method, but it captures the boundary (and I actually understand what it's doing LOL). 
Logged
billtavis
Safarist
******
Posts: 96


WWW
« Reply #14 on: February 20, 2015, 06:12:34 PM »

Here is a c++ function to get the gradient of a c-value based on the formulas given by Linas Vepstas. The gradient is returned by reference to the variables gzx and gzy. I also return the length and derivative values by reference to make it work with my other function below.
Code:
void getGradient(long double cx, long double cy, long double &xDer,
                 long double &yDer, long double &gzx, long double &gzy,
                 long double &length, unsigned int maxIter) {
        // based on gradient equation found at:
        // http://linas.org/art-gallery/escape/ray.html
    long double zx = 0.0L;
    long double zy = 0.0L;
    xDer = 0.0L; yDer = 0.0L; // initialize values
    long double xTemp;
    unsigned int n = 0;
    long double huge = std::numeric_limits<float>::max();
    while ((zx*zx + zy*zy) < huge && n < maxIter) {
            // compute new derivative:
        xTemp = 2.0L*(zx*xDer - zy*yDer) + 1.0L;
        yDer = 2.0L*(zx*yDer + zy*xDer);
        xDer = xTemp;
            // compute new z value:
        xTemp = zx*zx - zy*zy + cx;
        zy = 2.0L*zx*zy + cy;
        zx = xTemp;
        n++;
    }
    length = std::sqrt(zx*zx + zy*zy);
    long double m = n + 1 - (std::log(std::log(length)) / std::log(2.0L));
    long double f = pow(2.0L,-m); // Douady-Hubbard potential
    long double Dx = zx*(xDer/2.0) - zy*(-yDer/2.0); // real part of zn * Dn
    long double Dy = zx*(-yDer/2.0) + zy*(xDer/2.0); // imag part of zn * Dn
    gzx = f*Dx / (length*length*std::log(length)); // x-gradient for cx,cy
    gzy = f*Dy / (length*length*std::log(length)); // y-gradient for cx,cy
}

And here is a function that uses the getGradeint function to perform Runge-Kutta integration and find the external angle. It's assumed that the point has already been checked for being outside the set.
Code:
double externalAngle(long double cx, long double cy, std::size_t maxIter) {
        // uses 4th order Runge-Kutta integration to find external angle
    long double huge = 640000.0; // bigger is better, but also slower
    long double ckx = cx;
    long double cky = cy;
    long double ckxTemp, ckyTemp;
    long double xDer,yDer,length,derLength,dist,gzx,gzy,gzxA,gzyA,gzxB,gzyB,gzxC,gzyC,gzxD,gzyD,gzlength;
    unsigned int k = 0;
    unsigned int maxStep = maxIter*2; // make sure we try to walk far away
    while((ckx*ckx + cky*cky) < huge && k++ < maxStep) {
        // integrate with Runge-Kutta outwards along gradient:
        getGradient(ckx,cky,xDer,yDer,gzxA,gzyA,length,maxIter);
        derLength = std::sqrt(xDer*xDer + yDer*yDer);
        dist = 0.5L * std::log(length) * length / derLength;
        dist = std::min(16.0L,0.5*dist); // reduce step size
        gzlength = std::sqrt(gzxA*gzxA + gzyA*gzyA);
        ckxTemp = ckx + (0.5 * dist * gzxA/gzlength); // walk to midpoint
        ckyTemp = cky + (0.5 * dist * gzyA/gzlength); // using A
        getGradient(ckxTemp,ckyTemp,xDer,yDer,gzxB,gzyB,length,maxIter);
        gzlength = std::sqrt(gzxB*gzxB + gzyB*gzyB);
        ckxTemp = ckx + (0.5 * dist * gzxB/gzlength); // walk to midpoint
        ckyTemp = cky + (0.5 * dist * gzyB/gzlength); // using B
        getGradient(ckxTemp,ckyTemp,xDer,yDer,gzxC,gzyC,length,maxIter);
        gzlength = std::sqrt(gzxC*gzxC + gzyC*gzyC);
        ckxTemp = ckx + (dist * gzxC/gzlength); // walk full step
        ckyTemp = cky + (dist * gzyC/gzlength); // using C
        getGradient(ckxTemp,ckyTemp,xDer,yDer,gzxD,gzyD,length,maxIter);
        gzx = (1.0L/6.0L) * (gzxA + 2.0L*(gzxB + gzxC) + gzxD); // average the values together
        gzy = (1.0L/6.0L) * (gzyA + 2.0L*(gzyB + gzyC) + gzyD);
        gzlength = std::sqrt(gzx*gzx + gzy*gzy); // get length to normalize
        ckx += dist * gzx/gzlength; // normalize gradient and scale by dist
        cky += dist * gzy/gzlength;
    }
    return (std::atan2(cky,ckx));
}
Even though Runge-Kutta method requires solving the gradient 4 times per step, it winds up being much faster than Euler integration in practice because it is so much more accurate. Notice that I used
Code:
dist = std::min(16.0L,0.5*dist);
for the step size. To get an image just as smooth with plain Euler integration required reducing the step size to
Code:
dist = std::min(1.0L,0.05*dist);
which was way too slow.
« Last Edit: February 26, 2015, 05:29:12 PM by billtavis, Reason: fixed bug in code » Logged
Pages: [1] 2   Go Down
  Print  
 
Jump to:  


Powered by MySQL Powered by PHP Powered by SMF 1.1.21 | SMF © 2015, Simple Machines

Valid XHTML 1.0! Valid CSS! Dilber MC Theme by HarzeM
Page created in 0.142 seconds with 27 queries. (Pretty URLs adds 0.007s, 2q)