Can anyone help me understand interior distance estimation?

[laser blaster]

Hey. So, I've been trying to implement interior distance estimation into my fractal explorer, but I haven't been able to find a clear explanation on the interior DE algorithm.

Wikipedia has the formula on the Mandelbrot Set page, but the fancy notation used is completely indecipherable to me (I'm not a math expert, just a programmer).


The only code I've been able to find is here: http://www.moleculardensity.net/buddhabrot/appendix/2
Unfortunately, that code is so bloated and convoluted that I can't make any sense of it. It doesn't use many abstractions like vectors or functions, it's all written in terms of simple math operations on doubles. I just can't understand what it's doing.

Could someone with experience give me a simpler description of the interior DE algorithm? Thank you very much.

[Syntopia]

The interior distance estimate is quite complicated, but I actually think the link you provide gives quite a good description - formulas are shown both in component and complex notation.

If you have implemented the exterior estimate, some of the terms should be familiar. The biggest difference is that for interior estimate you have to keep track of four different running derivatives, instead of only one. In the code the four running derivaties are the D1,D2,D3,D4 variables.

You also have to find the periodic cycle, which may prove difficult if you are working in GLSL - I'm not sure you can allocate large enough float buffers.

Do you need the estimate? I don't think it will help you much with your AA issues, since there is no true link between the coverage and the distance estimate anyway.

[laser blaster]

Thank you for such a fast and thorough reply!

Hmm... I guess that code isn't really that difficult to understand. I was just intimidated by the bulk of it. The comments do explain what's going on in complex notation. I now understand that Z is the iterated point, and D1-D4 are the 4 derivatives mentioned on Wikipedia. The periodic cycle code is not as hard as I first thought, now that I realize it just returns a list of the points in the cycle that Z converges to.

But you're right, the periodic cycle code will be impractical to implement in glsl because it requires arbitrarily large buffers to store the points in the cycle. I'd say it's not very well-suited to GPU computing in general.

I don't exactly need the estimate, I was trying to use it to replace supersampling for AA (which can look decent, you just need to manually adjust a distance multiplier for every picture you take). But due to the potentially huge arrays required by the interior DE algorithm, I'm thinking it would be even more of a performance hog than SSAA.

Yeah, I guess I'll just have to stick with SSAA.

[therror]

Hello! Does anyone know if I can try it for the Multibrot set (with arbitrary non-integer powers)? I have implemented period detection for Multibrot and I want to add interior distance estimation for more complex coloring.

[Adam Majewski]

https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/demm#Interior_distance_estimation
hth

[therror]

Thank you very much, Adam! I have implemented Multibrot interior distance estimation! I think it is not fully correct. Maybe, period is computed with error.

[therror]

I fixed multiple errors, but unfortunately my formulae are wrong At power=2 all wright, but at higher powers (in this picture 7.97) only period=1 bulb (i.e. main part) shows correct distance estimation. Can anyone help me to find error?

Code:

  FLOAT zp, zq, dzp= ONE , dzq=0, dcp=0, dcq=0, dzdzp=0, dzdzq=0, dcdzp=0, dcdzq=0, dzpt, dzqt, dcpt, dcqt, dzdzpt, dcdzpt, x, y, cosat, sinat, x2y2=0, po, at, at2, sp, sq, xcompare, ycompare;\
....
  if (mindist<toler) mintest=2;\
  for (count=1; count<mintest; count++)\
  {\
    if (x2y2>8)\
    {\
      count=0;\
      break;\
    }\
    x2y2=x*x+y*y;\
    at2= ATAN2 (y, x);\
    po= POW (x2y2, HALF *powm);\
    at=powm*at2;\
    SINCOS (at, &sinat, &cosat);\
    x=po*cosat+sp;\
    y=po*sinat+sq;\
    po= POW (x2y2, HALF *(powm - ONE ));\
    at-=at2;\
    SINCOS (at, &sinat, &cosat);\
    zp=po*cosat;\
    zq=po*sinat;\
    dzpt=powm*(zp*dzp-zq*dzq);\
    dzqt=powm*(zq*dzp+zp*dzq);\
    dcpt=powm*(zp*dcp-zq*dcq)+ ONE ;\
    dcqt=powm*(zq*dcp+zp*dcq);\
    dzdzpt=powm*((zp*dzdzp-zq*dzdzq)+(dzp*dzp-dzq*dzq));\
    dzdzq =powm*((zp*dzdzq+zq*dzdzp)+(2*dzp*dzq));\
    dcdzpt=powm*((zp*dcdzp-zq*dcdzq)+(dzp*dcp-dzq*dcq));\
    dcdzq =powm*((zp*dcdzq+zq*dcdzp)+(dzp*dcq+dcp*dzq));\
    dzp=dzpt;\
    dzq=dzqt;\
    dcp=dcpt;\
    dcq=dcqt;\
    dzdzp=dzdzpt;\
    dcdzp=dcdzpt;\
    if ((xcompare-x)*(xcompare-x)+(ycompare-y)*(ycompare-y)<mindist) break;\
  }\

[therror]

Wrong second derivatives. Fixed.

Code:

for (count=1; count<mintest; count++)\
  {\
    if (x2y2>8)\
    {\
      count=0;\
      break;\
    }\
    x2y2=x*x+y*y;\
    at2= ATAN2 (y, x);\
    po= POW (x2y2, HALF *powm);\
    at=powm*at2;\
    SINCOS (at, &sinat, &cosat);\
    x=po*cosat+sp;\
    y=po*sinat+sq;\
    po= POW (x2y2, HALF *(powm - ONE ));\
    at-=at2;\
    SINCOS (at, &sinat, &cosat);\
    zp_1=po*cosat;\
    zq_1=po*sinat;\
    po= POW (x2y2, HALF *(powm - ONE - ONE ));\
    at-=at2;\
    SINCOS (at, &sinat, &cosat);\
    zp_2=po*cosat;\
    zq_2=po*sinat;\
    dzpt=powm*(zp_1*dzp-zq_1*dzq);\
    dzqt=powm*(zq_1*dzp+zp_1*dzq);\
    dcpt=powm*(zp_1*dcp-zq_1*dcq)+ ONE ;\
    dcqt=powm*(zq_1*dcp+zp_1*dcq);\
    dzp2=dzp*dzp-dzq*dzq;\
    dzq2=2*dzp*dzq;\
    dzdzpt=powm*((zp_1*dzdzp-zq_1*dzdzq)+(powm- ONE )*(zp_2*dzp2-zq_2*dzq2));\
    dzdzq =powm*((zp_1*dzdzq+zq_1*dzdzp)+(powm- ONE )*(zp_2*dzq2+zq_2*dzp2));\
    dzp2=dzp*dcp-dzq*dcq;\
    dzq2=dzp*dcq+dcp*dzq;\
    dcdzpt=powm*((zp_1*dcdzp-zq_1*dcdzq)+(powm- ONE )*(zp_2*dzp2-zq_2*dzq2));\
    dcdzq =powm*((zp_1*dcdzq+zq_1*dcdzp)+(powm- ONE )*(zp_2*dzq2+zq_2*dzp2));\
    dzp=dzpt;\
    dzq=dzqt;\
    dcp=dcpt;\
    dcq=dcqt;\
    dzdzp=dzdzpt;\
    dcdzp=dcdzpt;\
    if ((xcompare-x)*(xcompare-x)+(ycompare-y)*(ycompare-y)<mindist) break;\
  }\