### Riemandelettuce without trigonometry

• 33 Replies
• 858 Views

0 Members and 1 Guest are viewing this topic.

#### hobold #### Riemandelettuce without trigonometry

« 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);

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.

#### kosalos #### Riemandelettuce without trigonometry test

« 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..

#### hobold #### Re: Riemandelettuce without trigonometry

« 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/

#### mclarekin #### Re: Riemandelettuce without trigonometry

« 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

#### mclarekin #### Re: Riemandelettuce without trigonometry

« Reply #4 on: August 13, 2019, 01:47:21 PM »
power 8 is much better for rendering     #### mclarekin #### Re: Riemandelettuce without trigonometry

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

#### hobold #### Re: Riemandelettuce without trigonometry

« Reply #6 on: August 13, 2019, 02:13:54 PM »
power 8 is much better for rendering     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. #### claude #### Re: Riemandelettuce without trigonometry

« 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.

#### claude #### Re: Riemandelettuce without trigonometry

« 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).

#### claude #### Re: Riemandelettuce without trigonometry

« Reply #9 on: August 13, 2019, 03:43:45 PM »
another view

#### hobold #### Re: Riemandelettuce without trigonometry

« 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.

#### mclarekin #### Re: Riemandelettuce without trigonometry

« 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.

#### mclarekin #### Re: Riemandelettuce without trigonometry

« 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

#### hobold #### Re: Riemandelettuce without trigonometry

« 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;}

#### mclarekin #### Re: Riemandelettuce without trigonometry

« 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