Print Page - Can anyone help me understand interior distance estimation?

Title: Can anyone help me understand interior distance estimation?
Post by: laser blaster on April 04, 2013, 09:52:33 PM

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).
(http://upload.wikimedia.org/math/1/0/9/109e2d7ec490372864f5e67974ec5709.png)

The only code I've been able to find is here: http://www.moleculardensity.net/buddhabrot/appendix/2 (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.

Title: Re: Can anyone help me understand interior distance estimation?
Post by: Syntopia on April 04, 2013, 10:31:30 PM

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.

Title: Re: Can anyone help me understand interior distance estimation?
Post by: laser blaster on April 04, 2013, 11:43:54 PM

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.

Title: Re: Can anyone help me understand interior distance estimation?
Post by: therror on September 07, 2016, 09:32:24 AM

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.

Title: Re: Can anyone help me understand interior distance estimation?
Post by: Adam Majewski on September 07, 2016, 06:41:32 PM

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

Title: Re: Can anyone help me understand interior distance estimation?
Post by: therror on September 27, 2016, 11:42:28 AM

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.

Title: Re: Can anyone help me understand interior distance estimation?
Post by: therror on October 25, 2016, 03:15:33 PM

I fixed multiple errors, but unfortunately my formulae are wrong :sad1: 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;\
  }\

Title: Re: Can anyone help me understand interior distance estimation?
Post by: therror on October 29, 2016, 06:28:19 AM

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

Welcome to Fractal Forums

Fractal Software => Programming => Topic started by: laser blaster on April 04, 2013, 09:52:33 PM