Riemandelettuce without trigonometry

  • 33 Replies
  • 858 Views

0 Members and 1 Guest are viewing this topic.

Offline hobold

  • *
  • Fractal Friar
  • *
  • Posts: 115
« on: August 12, 2019, 08:48:20 AM »
I revisited my variant of the mandelbulb, the one based on monopolar spherical coordinates. I was trying to get a grasp of why the power 2 variant did not seem to exhibit an embedded 2D Mandelbrot Set in its relevant equator plane. (Despite the construction being based on this very embedding.)

Along the way, I replaced the trigonometric function calls of my older formulas with an application of the double angle formula. Computation is slightly less dog slow now. Here is a code fragment with the new computations:


// Riemandelettuce fractal iterations (monopolar spherical coordinates)
int iterateRmdltc(Vec3 p, int maxiter) {

  const double a = p.x;
  const double b = p.y;
  const double c = p.z;

  double len2 = p.dot(p);
  double rlen;
 
  int i = 0;
 
  while ((len2 < 4.0) && (i < maxiter)) {
   
    i++;

    rlen = 1.0/sqrt(len2);   
    p.scale(rlen);  // normalize vector to unit length => project onto sphere
   
    // find X-related iso-plane: polar projection onto unit circle
    double Kx = 2.0*p.x*(1.0 - p.y)/((p.y - 2.0)*p.y + p.x*p.x + 1.0);
    double Ky = 1.0 - 2.0*((p.y - 2.0)*p.y + 1.0) /
      ((p.y - 2.0)*p.y + p.x*p.x + 1.0);

    // doubled point
    double K2x = -2.0*Kx*Ky;
    double K2y = -(Ky*Ky - Kx*Kx);

    /*
    // two more doublings (for total power eight)
    Kx = -2.0*K2x*K2y;
    Ky = -(K2y*K2y - K2x*K2x);
    K2x = -2.0*Kx*Ky;
    K2y = -(Ky*Ky - Kx*Kx);
    */
   
    // (relevant) normal vector coordinates of doubled point plane
    double n1x = K2y - 1.0;
    double n1y = -K2x;

    // find Z-related iso-plane: polar projection onto unit circle
    double Kz = 2.0*p.z*(1.0 - p.y)/((p.y - 2.0)*p.y + p.z*p.z + 1.0);
    Ky = 1.0 - 2.0*((p.y - 2.0)*p.y + 1.0)/((p.y - 2.0)*p.y + p.z*p.z + 1.0);

    // doubled point
    double K2z = -2.0*Kz*Ky;
    K2y = -(Ky*Ky - Kz*Kz);

    /*
    // two more doublings (for total power eight)
    Kz = -2.0*K2z*K2y;
    Ky = -(K2y*K2y - K2z*K2z);
    K2z = -2.0*Kz*Ky;
    K2y = -(Ky*Ky - Kz*Kz);
    */
   
    // (relevant) normal vector coordinates of doubled point plane
    double n2y = -K2z;
    double n2z = K2y - 1.0;

    // compute position of doubled point as intersection of planes and sphere
    // solved ray parameter
    double nt = 2.0*(n1x*n1x*n2z*n2z)/((n1x*n1x + n1y*n1y)*n2z*n2z
                   + n1x*n1x*n2y*n2y);
       
    // doubled point position
    p.y = 1.0 - nt;
    p.x = n1y*(1.0 - p.y)/n1x;
    p.z = n2y*(1.0 - p.y)/n2z;
   
    // raise original length to the power, then add constant
    //p.scale(len2*len2*len2*len2);  // for 8th power
    p.scale(len2);
    p.add(Vec3(a,b,c));
   
    len2 = p.dot(p);
  }
   
  if (len2 < 4.0) {
    return -1;
  }

  return i;
}


I still don't have a closed formula. I strongly suspect one would end up with a rational function that has quartic polynomials in both enumerator and denominator. But that means up to 35 coefficients in either polynomial, and the symbolic formalism quickly got out of hand.

The elusiveness of the embedded 2D Mandelbrot Set seems to be partly caused by numerical instability; there is actually a division by zero in all of the "equator plane". I am currently not sure if that was just a bad choice of parameterization on my part, or if the denominator really has to be zero in that spot. Geometrically, though, the embedded plane does not behave very differently than nearby points; there should be no real discontinuity there.

Linkback: https://fractalforums.org/fractal-mathematics-and-new-theories/28/riemandelettuce-without-trigonometry/2996/

Offline kosalos

  • *
  • Fractal Friar
  • *
  • Posts: 102
« Reply #1 on: August 12, 2019, 06:52:42 PM »
I used your equation as my distance estimation function,
with 'two more doublings' enabled, and  log(length(p))  as DE value (iteration count unused).
My lettuce doesn't do too much.
Like to see a snapshot of what it SHOULD look like..

Offline hobold

  • *
  • Fractal Friar
  • *
  • Posts: 115
« Reply #2 on: August 12, 2019, 10:27:38 PM »
The presented formula is not a distance estimate. It only computes the orbit of a point. So the result means "inside" (not escaped) when -1 is returned, otherwise the number of iterations until bailout.

The power 8 set should look somewhat like this (top post on page):
http://www.fractalforums.com/the-3d-mandelbulb/riemandelettuce/30/

Offline mclarekin

  • *
  • Fractal Frankfurter
  • *
  • Posts: 612
« Reply #3 on: August 13, 2019, 01:34:53 PM »
a quick attempt, without trying analytic DE (which from memory is tricky ) top image is a julia

Offline mclarekin

  • *
  • Fractal Frankfurter
  • *
  • Posts: 612
« Reply #4 on: August 13, 2019, 01:47:21 PM »
power 8 is much better for rendering  :) :) O0 O0 O0

Offline mclarekin

  • *
  • Fractal Frankfurter
  • *
  • Posts: 612
« Reply #5 on: August 13, 2019, 02:00:34 PM »
y axis ( these have a scale of 0.7 applied)

Offline hobold

  • *
  • Fractal Friar
  • *
  • Posts: 115
« Reply #6 on: August 13, 2019, 02:13:54 PM »
power 8 is much better for rendering  :) :) O0 O0 O0
So you managed to pull a distance estimator out of the hat given the new formulation? I should have done this work years ago ...

And yes, the power 2 Riemandelettuce is very "leafy", there are many thin and long sheets protruding in many directions. It's a lettuce after all. :)

Offline claude

  • *
  • 3f
  • ******
  • Posts: 1273
    • mathr.co.uk
« Reply #7 on: August 13, 2019, 02:38:55 PM »
Distance estimate is related to inverse of magnitude of gradient of continuous escape count field. This can be calculated numerically (eg by sampling at several nearby points) rather than by using analytic derivatives.

More precisely, for quadratic Mandelbrot set, if estimate d = 1/(log(2) |dn/dc|) is small enough then true distance is between d/(2+e) .. (2+e)*d for small e.

Offline claude

  • *
  • 3f
  • ******
  • Posts: 1273
    • mathr.co.uk
« Reply #8 on: August 13, 2019, 03:32:03 PM »
here's an attempt at porting the power 8 code to FragM using Raymond:

Code: [Select]
float Lettuce(vec3 p, int maxiter)
{
  float a = p.x;
  float b = p.y;
  float c = p.z;
  float len2 = dot(p, p);
  float rlen;
  int i = 0;
  while ((len2 < 4.0) && (i < maxiter)) {
    i++;
    rlen = 1.0/sqrt(len2);   
    p *= rlen;  // normalize vector to unit length => project onto sphere
    // find X-related iso-plane: polar projection onto unit circle
    float Kx = 2.0*p.x*(1.0 - p.y)/((p.y - 2.0)*p.y + p.x*p.x + 1.0);
    float Ky = 1.0 - 2.0*((p.y - 2.0)*p.y + 1.0) /
      ((p.y - 2.0)*p.y + p.x*p.x + 1.0);
    // doubled point
    float K2x = -2.0*Kx*Ky;
    float K2y = -(Ky*Ky - Kx*Kx);
    // two more doublings (for total power eight)
    Kx = -2.0*K2x*K2y;
    Ky = -(K2y*K2y - K2x*K2x);
    K2x = -2.0*Kx*Ky;
    K2y = -(Ky*Ky - Kx*Kx);
    // (relevant) normal vector coordinates of doubled point plane
    float n1x = K2y - 1.0;
    float n1y = -K2x;
    // find Z-related iso-plane: polar projection onto unit circle
    float Kz = 2.0*p.z*(1.0 - p.y)/((p.y - 2.0)*p.y + p.z*p.z + 1.0);
    Ky = 1.0 - 2.0*((p.y - 2.0)*p.y + 1.0)/((p.y - 2.0)*p.y + p.z*p.z + 1.0);
    // doubled point
    float K2z = -2.0*Kz*Ky;
    K2y = -(Ky*Ky - Kz*Kz);
    // two more doublings (for total power eight)
    Kz = -2.0*K2z*K2y;
    Ky = -(K2y*K2y - K2z*K2z);
    K2z = -2.0*Kz*Ky;
    K2y = -(Ky*Ky - Kz*Kz);
    // (relevant) normal vector coordinates of doubled point plane
    float n2y = -K2z;
    float n2z = K2y - 1.0;
    // compute position of doubled point as intersection of planes and sphere
    // solved ray parameter
    float nt = 2.0*(n1x*n1x*n2z*n2z)/((n1x*n1x + n1y*n1y)*n2z*n2z
                   + n1x*n1x*n2y*n2y);
    // doubled point position
    p.y = 1.0 - nt;
    p.x = n1y*(1.0 - p.y)/n1x;
    p.z = n2y*(1.0 - p.y)/n2z;
    // raise original length to the power, then add constant
    // p *= len2;
    p *= len2*len2*len2*len2;  // for 8th power
    p += vec3(a,b,c);
    len2 = dot(p, p);
  }
  if (len2 < 4.0) {
    return -1.0;
  }
  return float(i) + 1 - log(sqrt(len2)) / log(8.0);
}

Then using numerically calculated spatial gradient over nearby points for de and surface normal (should really use gradient of de for that but it's slow enough as-is).

Offline claude

  • *
  • 3f
  • ******
  • Posts: 1273
    • mathr.co.uk
« Reply #9 on: August 13, 2019, 03:43:45 PM »
another view

Offline hobold

  • *
  • Fractal Friar
  • *
  • Posts: 115
« Reply #10 on: August 13, 2019, 09:49:25 PM »
By the way, that image of a julia set above looks quite ... harmonic. I don't know, like, crystallized music? Certainly not fractally rough. Quite unexpected.

Offline mclarekin

  • *
  • Fractal Frankfurter
  • *
  • Posts: 612
« Reply #11 on: August 14, 2019, 05:10:51 AM »
@ hobold, i think it is because the julia is clipped by the bailout radius??

@claude. I used Buddhi's deltaDE for these images, which i think is what you saying. I find with these riemann spheres i have to include a scale below 1.0 to get the DeltaDE to work well.  (Normally 0.7 to 0.95)

I can get this one to run with an analytic fudge,  I will see if i can fine tune it.

Offline mclarekin

  • *
  • Fractal Frankfurter
  • *
  • Posts: 612
« Reply #12 on: August 15, 2019, 01:43:57 PM »
julias
second row some hybrids.

I got some form of analytic DE working but  detail gets lost

Offline hobold

  • *
  • Fractal Friar
  • *
  • Posts: 115
« Reply #13 on: August 20, 2019, 10:41:11 AM »
The word "hybrid" and an unrelated question about DE formulas inspired me to try another experiment: hybrid exponent lettuce. Have a look (82MByte MP4 download):

http://vectorizer.org/rmdltc/pow6_4_4_iter10.mp4

The monopolar spherical coordinates do not have distinguishable longitude and latitude. Instead, the two coordinate components are interchangeable and essentially the same. One can use different wraparound exponents for either!

The above animation loop uses 6 and 4 as wraparound exponents, and 4 again as stretch exponent. IMHO, the variety and diversity of shapes is notably richer as compared to using the same wraparound exponent on both coordinates.

Here is the code again for complete specification:


int iterateRmdltc(Vec3 p, int maxiter) {

  const double a = p.x;
  const double b = p.y;
  const double c = p.z; 

  double len2 = p.dot(p);
  double rlen;
 
  int i = 0;
 
  while ((len2 < 4.0) && (i < maxiter)) {
   
    i++;

    rlen = 1.0/sqrt(len2);   
    p.scale(rlen);  // normalize vector to unit length => project onto sphere
   
    // find X-related iso-plane: polar projection onto unit circle
    double Kx = 2.0*p.x*(1.0 - p.y)/((p.y - 2.0)*p.y + p.x*p.x + 1.0);
    double Ky = 1.0 - 2.0*((p.y - 2.0)*p.y + 1.0) /
      ((p.y - 2.0)*p.y + p.x*p.x + 1.0);

    // doubled point
    double K2x = -2.0*Kx*Ky;
    double K2y = -(Ky*Ky - Kx*Kx);
   
    // one more doubling (total power 4)
    Kx = -2.0*K2x*K2y;
    Ky = -(K2y*K2y - K2x*K2x);
    K2x = Kx;
    K2y = Ky;   
   
    // (relevant) normal vector coordinates of doubled point plane
    double n1x = K2y - 1.0;
    double n1y = -K2x;

    // find Z-related iso-plane: polar projection onto unit circle
    double Kz = 2.0*p.z*(1.0 - p.y)/((p.y - 2.0)*p.y + p.z*p.z + 1.0);
    Ky = 1.0 - 2.0*((p.y - 2.0)*p.y + 1.0)/((p.y - 2.0)*p.y + p.z*p.z + 1.0);

    // doubled point
    double K2z = -2.0*Kz*Ky;
    K2y = -(Ky*Ky - Kz*Kz);

    // another tripling (total power 6)
    Kz = -K2z*(3.0 - 4.0*K2z*K2z);
    Ky = -K2y*(4.0*K2y*K2y - 3.0);
    K2x = Kx;
    K2y = Ky;

    // (relevant) normal vector coordinates of doubled point plane
    double n2y = -K2z;
    double n2z = K2y - 1.0;

    // compute position of doubled point as intersection of planes and sphere
    // solved ray parameter
    double nt = 2.0*(n1x*n1x*n2z*n2z)/((n1x*n1x + n1y*n1y)*n2z*n2z
                   + n1x*n1x*n2y*n2y);
       
    // warped point position
    p.y = 1.0 - nt;
    p.x = n1y*(1.0 - p.y)/n1x;
    p.z = n2y*(1.0 - p.y)/n2z;
   
    // raise original length to the power, then add constant
    p.scale(len2*len2);  // 4th power
    p.add(Vec3(a,b,c));
   
    len2 = p.dot(p);
  }
   
  if (len2 < 4.0) {
    return -1;
  }

  return i;
}


Offline mclarekin

  • *
  • Fractal Frankfurter
  • *
  • Posts: 612
« Reply #14 on: August 20, 2019, 01:01:46 PM »
cool, thanks, i will try and implement, plus the pow4

attached is another pow8 hybrid image