(Question) analytic normals for inflated Mandelbrot set?

  • 14 Replies
  • 286 Views

0 Members and 1 Guest are viewing this topic.

Offline claude

  • *
  • 3c
  • ***
  • Posts: 895
    • mathr.co.uk
« on: November 12, 2018, 08:45:30 PM »
Consider the interior distance estimate for the Mandelbrot set:
https://en.wikipedia.org/wiki/Mandelbrot_set#Interior_distance_estimation
Using the interior distance as height gives cone-like shapes.  Modifying the numerator to
\[ \sqrt{1 - \left|\frac{\partial}{\partial z}\right|^2} \]
gives more sphere-like shapes, this works as the numerator is in [0..1] for interior points.

The problem now is, given a point C in the interior of a component of period P, with the surface height as described above, how to calculate the surface normal in the object space relative to C.

One tangent in the component space (conformally mapped to the unit disc by \( \frac{\partial}{\partial z} \)) is in the direction \( \left(i \frac{\partial}{\partial z}, 0\right) \) but I'm not sure how to unmap that to be relative to C.  The other tangent is harder to find, and both are needed for the vector cross product to find the normal (unless there is a more direct way?).

The answer would be nice, but explanation would help too, as I would like to do the same thing for atom domains later.

Offline FractalDave

  • *
  • Uploader
  • *
  • Posts: 167
    • Makin Magic Fractals
« Reply #1 on: November 12, 2018, 11:29:27 PM »
For outside DE angles I stuck to fully calculating in complex all the way through and used the resultant angle of the complex DE value, though I only used it as a colouring not actually for normals and I never tried either a full W/N triplex or a bicomplex or a quaternion version for higher dimensions = the colouring does look promising though, at least for "outside" normals, I don't know if a similar method could be applied for "inside" DE or not.
The attached image of the DE angles on the z^2+c Mandelbrot was created in Ultra Fractal, parameters below, all formulas are in the UF public formula database, requires at least UF5 I think.

Code: [Select]
Fractal1 {
::EK3KKjn2deVTTuJOQ07uK/fgi7xGhHw4kSH2NTSusTukc3lGJhHlRIRJEOj3f9pBBe4DRGnc
  xlRv+186+JJkyNEqlIf/6VBBWhVyxhf2NEKM4nCm9JMK9QawTcxpns4DJ3FIJX4mKMK4Epog
  gj3k0Qma4MhtCHeP5sgF8A5Zh6DI0WU824IUWYARZFEpgUJUnwX4VrX1mo23MlUaFaFO8fJ0
  nPZ01KWYgukQF2LYUU06VFkySgoLauyyN4oNxZ7PsHFjQgs2v9dRbSju7QcaW6hd73hSCKIn
  U49JI0mskMQl5aTRtk0mkCyLimsAZ/QQuQyVkCo4/mloYEDbTdeRYA8iMXwhPAjxlPa02wgy
  jVWixij2G1kmyjl6fCpJebE8/HJCpu2iTTS2luelQVJYcXv1QUV5QgKtiveFE0cIpQxJmBip
  oIPZTNVG2EWnW+CXxNC63K0a7TfxQYCA4jap2AtnW51i00kvXcmbOBwNj+soEH7U8IVCPzI1
  WdjlMpayrZn4Y0moO+UMylgqCiU2kFM6TvLBGgOID0Zpg2lDau4FOD6Aus4e8MBySfvk2VHf
  USqqar/N1yHf/DP85j9Gz9Co/ro8PVZFFEr2E6h6mzHP56TN+rH0+ndTu8yvr1+2Ry6UEvXR
  PSq4LHOUUNThaLmZUhOFyHth90xI/oGWT9aPdML4HOzepkDBsAapWos/W+l6KR7qTP5YyE+5
  pgpZiz9zMGhUoZc8/oOJ5e5JOXplCmXqAIsIRi92hHN5Mzjgoa15FKYaro6niNHGYeVVe5DB
  4U28eJ3cVU+bUlt7OMn4IN5jKhxAuCNzH9GwXNQvuM30uFEOkqLK0q2lcftWK/e3GTh/OWjX
  q9GRedru3OUFIgbIcoq/7EfPxbS/DC+PQT3aVor+bKhOW3g+fNybVK3e/f5FpA41lCwkw2P2
  2PBsPxe2nHK6yJF+oont1+U05ieaEN7ADjVOMyZhVeUnnXxtRYvZp84/XRJSeEcOoloiWiKc
  OEKH61u1rTRfpRS/X7ZB8TXQsdfPwDoD5O/Md6KGHtY9g6/Gvfq76LpG7s3meTDcYwz8vJgz
  tvJBsg7NKqZm3E0henXioFIO25mAOz4mgPy3mgNy28KJnr5tQuaaeJe1zyFKi8mNtRRPz1mi
  O32mGxC+24wmZcThH6c+piWi6YvbK6MzbaAjcvpgjsP/6y5f+rnrGofqtOY/BSbvzg7Q+t77
  Kg7m8Ck72dfxZ7RpxZZdj20idjjSh7Jllk0BE7cJH2dpJHQdA7elRCcguu7gN8dCdv+XZ/N0
  iTgLX9LAA3OB0B==
}

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

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1536
« Reply #2 on: November 13, 2018, 12:17:06 AM »
http://mathworld.wolfram.com/NormalVector.html does not answer your question?
Maybe I don't get "normal in object space relative to c". Something else than the surface normal?

Offline claude

  • *
  • 3c
  • ***
  • Posts: 895
    • mathr.co.uk
« Reply #3 on: November 13, 2018, 05:56:29 PM »
Thanks gerrit, it helps a bit, but the expression for \( f(x, y) \) is quite hairy, with part of it being implicit:

Given complex \( C \) in an atom domain of period \( P \), find \( Z_0 = F_C^P(Z_0) \) where \( F_c(z)=z^2+c \).  Now, if \( \left|\frac{\partial}{\partial z}F_C^P(Z_0)\right| \le 1 \) then \( C \) is interior to a component of period \( P \).  Now define the surface height
\[ h = f\left(\Re(C), \Im(C)\right) = \frac{ \sqrt{1 - \left| \frac{\partial}{\partial z} \right|^2} }{\left|\frac{\partial}{\partial{c}}\frac{\partial}{\partial{z}} + \frac{\left(\frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}\right) \left(\frac{\partial}{\partial{c}}\right)} {1-\frac{\partial}{\partial{z}}}\right|} \]
with all the derivatives in that expression being evaluated at \( F_C^P(Z_0) \).

Finding the derivatives of \( f \) seems to be a big challenge, I was wondering if maybe there is a simpler way using properties of conformal mapping etc?

Offline claude

  • *
  • 3c
  • ***
  • Posts: 895
    • mathr.co.uk
« Reply #4 on: November 13, 2018, 06:06:48 PM »
btw this is my work in progress with bad normals (simply the direction to the nucleus of the component)

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1536
« Reply #5 on: November 14, 2018, 01:39:40 AM »
Thanks gerrit, it helps a bit, but the expression for \( f(x, y) \) is quite hairy, with part of it being implicit:

Given complex \( C \) in an atom domain of period \( P \), find \( Z_0 = F_C^P(Z_0) \) where \( F_c(z)=z^2+c \).  Now, if \( \left|\frac{\partial}{\partial z}F_C^P(Z_0)\right| \le 1 \) then \( C \) is interior to a component of period \( P \).  Now define the surface height
\[ h = f\left(\Re(C), \Im(C)\right) = \frac{ \sqrt{1 - \left| \frac{\partial}{\partial z} \right|^2} }{\left|\frac{\partial}{\partial{c}}\frac{\partial}{\partial{z}} + \frac{\left(\frac{\partial}{\partial{z}}\frac{\partial}{\partial{z}}\right) \left(\frac{\partial}{\partial{c}}\right)} {1-\frac{\partial}{\partial{z}}}\right|} \]
with all the derivatives in that expression being evaluated at \( F_C^P(Z_0) \).

Finding the derivatives of \( f \) seems to be a big challenge, I was wondering if maybe there is a simpler way using properties of conformal mapping etc?
OK, note that h^2 is the sum of two absolute values of analytic functions, and of course \( \nabla h = \nabla h^2 /(2h) \).
Section 2 of http://persianney.com/fractal/fractalNotes.pdf shows how to go back and forth between \( R^2 \) and \( C \), both for variables and gradient/derivatives.

For analytic function f() consider grad of |f|. In the notes I show that norm of grad is just abs of derivative of f. Similarly you can work out that
\( (\partial_x - i \partial_y)|f| = \partial_z|f| = \bar{f}f^{'}/(2|f|) \). So all you have to work out is \( \partial_c  h^2 \) which is only a bit unpleasant.

Offline claude

  • *
  • 3c
  • ***
  • Posts: 895
    • mathr.co.uk
« Reply #6 on: November 14, 2018, 05:19:34 PM »
Thanks that helps a lot! I get:
\[ h^2 = |A| - |B| = \left|\frac{1}{a}\right| - \left|\frac{b}{a}\right| \]
\[ \frac{\partial}{\partial c} h^2 = \frac{\overline{A} \frac{\partial}{\partial c} A}{2 |A|} - \frac{\overline{B} \frac{\partial}{\partial c} B}{2 |B|} \]
\[ \frac{\partial}{\partial c} A = -\frac{\frac{\partial}{\partial c} a}{a^2} \]
\[ \frac{\partial}{\partial c} B = \frac{\frac{\partial}{\partial c} b}{a}-\frac{\frac{\partial}{\partial c} a}{a^2} \]
which all means that \( \frac{\partial}{\partial c} h^2 \) can be found without too much further trouble (an A4 page of equations in my not so small handwriting results from wolfram alpha, probably doing it in stages in numerical code would make more sense than expanding everything, anyway)
plug that and h into your \(  \nabla h = \nabla h^2 / 2h \) and that's the result, using \( \nabla F = \left(\Re{\frac{\partial}{\partial c} F}, -\Im{\frac{\partial}{\partial c} F}\right) \)

Offline claude

  • *
  • 3c
  • ***
  • Posts: 895
    • mathr.co.uk
« Reply #7 on: November 15, 2018, 01:12:05 PM »
something is still not right - maybe I made a typo?

Code: [Select]
struct Mandelbrot_interior
{
  vec2 dz;
  vec3 N;
  float y;
};

Mandelbrot_interior Mandelbrot_interior_de(vec2 c, int period, vec2 z0)
{
        vec2 w = z0;
        vec2 dz = vec2(1.0, 0.0);
        vec2 dc = vec2(0.0);
        vec2 dzdz = vec2(0.0);
        vec2 dcdz = vec2(0.0);
        vec2 dcdc = vec2(0.0);
        vec2 dcdzdz = vec2(0.0);
        vec2 dcdcdz = vec2(0.0);
        for (int j = 0; j < period; ++j)
        {
          dcdcdz = 2.0 * (cMul(w, dcdcdz) + 2.0 * cMul(dc, dcdz) + cMul(dz, dcdc));
          dcdzdz = 2.0 * (cMul(w, dcdzdz) + cMul(dc, dzdz) + 2.0 * cMul(dcdz, dz));
          dcdc = 2.0 * (cMul(w, dcdc) + cSqr(dc));
          dcdz = 2.0 * (cMul(w, dcdz) + cMul(dz, dc));
          dzdz = 2.0 * (cMul(w, dzdz) + cSqr(dz));
          dc = 2.0 * cMul(w, dc) + vec2(1.0, 0.0);
          dz = 2.0 * cMul(w, dz);
          w = cSqr(w) + c;
        }
        vec2 a = dcdz + cDiv(cMul(dzdz, dc), vec2(1.0,0.0) - dz);
        float h = sqrt(1.0 - dot(dz, dz)) / length(a);
        vec2 Da = dcdcdz + cDiv(cMul(dcdzdz, dc) + cMul(dzdz, dcdc), vec2(1.0,0.0) - dz)
                         + cDiv(cMul(cMul(dzdz, dc), dcdz), cSqr(vec2(1.0,0.0) - dz));
        vec2 dA = -cDiv(Da, cSqr(a));
        vec2 A = cDiv(vec2(1.0, 0.0), a);
        vec2 b = cSqr(dz);
        vec2 Db = 2.0 * cMul(dz, dcdz);
        vec2 dB = -cDiv(cMul(b, Da) - cMul(Db, a), cSqr(a));
        vec2 B = cDiv(b, a);
        vec2 dh2 = cMul(cConj(A), dA) / (2.0 * length(A))
                 - cMul(cConj(B), dB) / (2.0 * length(B));
        vec2 gradh = cConj(dh2 / (2.0 * h));
        vec3 N = vec3(gradh, -1.0);
        Mandelbrot_interior r;
        r.dz = dz;
        r.N = N;
        r.y = h;
        return r;
}

...
      Mandelbrot_interior d = Mandelbrot_interior_de(P.xy, period, Mandelbrot_attractor(P.xy, period));
      vec3 N = normalize(d.N);
      if (dot(N, direction) < 0.0) N = -N;
...

attached image is with diffuse surface material, the bulbs look bad...

Offline claude

  • *
  • 3c
  • ***
  • Posts: 895
    • mathr.co.uk
« Reply #8 on: November 15, 2018, 06:35:32 PM »
nevermind about the normals, my algorithm is flawed at a more basic level in 2 ways:

1. (fixable) when checking that points are exterior to a component, need to check lower periods (dividing the target period) too in case it is interior to one of instead..
2. (fatal?) the nearest point on the surface of a component may not be on the straight line to the nucleus (especially for cardioids), so finding the distance to this point as a DE is meaningless

(1) causes the "collar" around the bulb connection (afaict)
(2) causes farting from the cardioid cusp, and probably other badness...

I go back to thinking hard about this problem...

Offline knighty

  • *
  • Fractal Feline
  • **
  • Posts: 185
« Reply #9 on: November 15, 2018, 06:48:44 PM »
Looks cool! :)

Offline claude

  • *
  • 3c
  • ***
  • Posts: 895
    • mathr.co.uk
« Reply #10 on: November 16, 2018, 04:03:05 AM »
I improved my nucleus-finding algorithm, to reject lower periods automatically, thus enlarging the Newton basins and making it much more likely to converge to the right place (added lines marked /* */ in the below which uses Complex.frag dual numbers):

Code: [Select]
vec2 Mandelbrot_nucleus(vec2 c0, int period)
{
        vec2 c = c0;
        // 8 Newton iterations is too few, 16 seems better
        for (int j = 0; j < 16; ++j)
        {
          vec4 G = cConst(1.0);
          vec4 z = cConst(0.0);
          for (int l = 1; l <= period; ++l)
          {
            z = cSqr(z) + cVar(c);
/* */       if (l < period && period % l == 0)
/* */         G = cMul(z, G);
          }
          G = cDiv(z, G);
          c -= cDiv(G.xy, G.zw);
        }
        return c;
}

I improved my attractor implementation in the same way (it returns the derivative as well as z_0 because it's useful for other algorithms like interior checking):

Code: [Select]
vec4 Mandelbrot_attractor(vec2 c, vec2 z_initial, int period)
{
        vec4 z = cVar(z_initial);
        vec2 dz = vec2(0.0);
        for (int j = 0; j < 16; ++j)
        {
          vec4 z0 = z;
          vec4 G = cConst(1.0);
          for (int l = 1; l <= period; ++l)
          {
            z = cSqr(z) + cConst(c);
/* */       if (l < period && period % l == 0)
/* */         G = cMul(z - z0, G);
          }
          dz = z.zw;
          G = cDiv(z - z0, G);
          z = cVar(z0.xy - cDiv(G.xy, G.zw));
        }
        return vec4(z.xy, dz);
}

I also now trace the attractor along a ray from the nucleus, which improves stability in the exterior:

Code: [Select]
vec4 Mandelbrot_attractor_ray(vec2 nucleus, vec2 target, int period)
{
  const int steps = 64;
  vec4 P = vec4(0.0);
  for (int ray = 1; ray <= steps; ++ray)
  {
    P = Mandelbrot_attractor(mix(nucleus, target, float(ray) / float(steps)), P.xy, period);
  }
  return P;
}

Attached are comparisons of the period 4 basins with the two methods (naive unmodified attractor (20) vs modified everything with attractor_ray (19)).

Offline claude

  • *
  • 3c
  • ***
  • Posts: 895
    • mathr.co.uk
« Reply #11 on: November 16, 2018, 04:14:20 AM »
for sphere-like inflated disc-like shapes, this works ok:

Code: [Select]
...
      c = pos.xy;
      vec2 direction = normalize(dz_of_c);
      vec2 delta = c - nucleus;
      float tanthetasquared = (pos.z * pos.z) / dot(delta, delta);
      float probe = 1.0 / sqrt(1.0 + tanthetasquared);
      vec2 c_target = Mandelbrot_interior_point(nucleus, probe * direction, period);
...

where probe is the cosine of the angle from horizontal at the nucleus of the input position, indicating how far out from the nucleus the intersection is (ie, the surface height (relative to component radius) is the sine of this angle, assuming perfect sphere).  should do some refinement passes to find the closest point, but it seems ok enough for now.

cardioids are still woeful.  reducing fudge factor makes it more solid, but expands everything else too - looks bad too.

Offline mclarekin

  • *
  • Fractal Fluff
  • *****
  • Posts: 378
« Reply #12 on: November 16, 2018, 07:35:53 AM »
top image is really beautiful O0
bottom is a good example of fudge factor fattening :)

Offline claude

  • *
  • 3c
  • ***
  • Posts: 895
    • mathr.co.uk
« Reply #13 on: November 16, 2018, 03:32:53 PM »
I think I solved it by scrapping the interior distance estimate idea and using something more calculable instead:
Code: [Select]
...
      c  = object.xy;
      nucleus = Mandelbrot_nucleus(c, period);
      if (period > 1 && rand < 0.2)      {
        atom = cDiv(z, mz);
        S = Mandelbrot_domain_size(nucleus, period);
      }      else      {
        atom = Mandelbrot_attractor(nucleus, c, period).dz;
        S = 0.5 * length(Mandelbrot_size(nucleus, period));
      }
      if (S > 0.001)      {
        d = (length(vec3(atom, object.z / S)) - 1.0) * S;
        if (d < de) de = d;
      }
      if (length(z) < length(mz)) mz = z;
...

added transparent atom domains by discarding them stochastically.

but slow shader is slow, even with the default simple DE-Raytracer.frag

Offline knighty

  • *
  • Fractal Feline
  • **
  • Posts: 185
« Reply #14 on: November 16, 2018, 06:23:32 PM »
 :thumbs:
I wasn't expecting it to work so well on GPU because of the limited FP precision.


xx
"Time Span"

Started by cricke49 on Fractal Image Gallery

0 Replies
119 Views
Last post August 02, 2018, 07:05:21 AM
by cricke49
xx
Birdie Style

Started by gannjondal on Fractal Image Gallery

1 Replies
210 Views
Last post May 08, 2018, 02:39:37 PM
by who8mypnuts
xx
Buddhabrot-style Burning Ship [65536x24576]

Started by programagor on Fractal Image Gallery

12 Replies
260 Views
Last post October 10, 2018, 02:12:39 AM
by 3DickUlus
clip
analytic DE for skewed images

Started by claude on Kalles Fraktaler

4 Replies
175 Views
Last post April 19, 2018, 08:50:39 AM
by claude
clip
analytic log(DE) post-processed with ImageMagick

Started by claude on Fractal Image Gallery

3 Replies
129 Views
Last post April 03, 2018, 07:00:04 AM
by claude