Logo by mario837 - Contribute your own Logo!
News: Support us via Flattr FLATTR Link
 
*
Welcome, Guest. Please login or register. September 25, 2017, 05:17:55 PM


Login with username, password and session length


The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!


Pages: 1 2 [3]   Go Down
  Print  
Share this topic on DiggShare this topic on FacebookShare this topic on GoogleShare this topic on RedditShare this topic on StumbleUponShare this topic on Twitter
Author Topic: A new? Simple way to compute DE for any trig. Mandelbulb  (Read 7064 times)
0 Members and 1 Guest are viewing this topic.
Buddhi
Fractal Iambus
***
Posts: 885



WWW
« Reply #30 on: January 09, 2010, 04:53:41 PM »

double th = atan2(y, x);
double ph = -atan2(z, sqrt(x * x + y * y));
x = r2 * cos(p * ph) * cos(p * th) + par.a;
y = r2 * cos(p * ph) * sin(p * th) + par.b;
z = r2 * sin(p * ph) + par.c;
Logged

David Makin
Global Moderator
Fractal Senior
******
Posts: 2286



Makin' Magic Fractals
WWW
« Reply #31 on: January 09, 2010, 05:46:36 PM »

How about timing? This 640*640 image took 35 seconds to render onscreen (so no anti-aliasing)
Using my program it takes 20 seconds to render 640x640 image (step distance was 0.3*DE).

Hi guys, just to point out that if you're comparing timings you really need to state the computing power used wink


Managed to get a similar view of the same Julia...



Rendered within UF @640*640 in 27secs on my 3GHz P4HT using the current optimised version of my formula, just to put that in context exactly the same image rendered using full trig for both the iterate and the derivative takes nearly 5 minutes.
Rendered using analytical DE, threshold 1e-3, step sizes = 0.45*DE.

Here's part of the optimised code:

              magn = sqrt(magn)
              dr = @mpwr*(ph=sqrt((th=|dzri|) + sqr(dzj)))
              if (r=|zri|)>0.0
                zjk = ((r=sqrt(r)) + flip(zj))/magn
                dr = dr*magn^(@mpwr-1.0)
                zri = zri/r
                if @mpwr==round(@mpwr) && abs(@mpwr-1.0)<32
                  if @mpwr==1.0
                    zri1 = (1,0)
                    zjk1 = (1,0)
                  elseif @mpwr==0.0
                    zri1 = 1.0/zri
                    zjk1 = 1.0/zjk
                  else
                    if @mpwr-1.0<0.0
                      zri1 = 1.0/zri
                      zjk1 = 1.0/zjk
                    else
                      zri1 = zri
                      zjk1 = zjk
                    endif
                    if abs(@mpwr-1.0)==2.0
                      zri1 = sqr(zri1)
                      zjk1 = sqr(zjk1)
                    elseif abs(@mpwr-1.0)==3.0
                      zri1 = zri1^3
                      zjk1 = zjk1^3
                    elseif abs(@mpwr-1.0)==4.0
                      zri1 = zri1^4
                      zjk1 = zjk1^4
                    else
                      if (abs(@mpwr-1.0)%4)>1
                        zp2 = sqr(zri1)
                        z1p2 = sqr(zjk1)
                      endif
                      if (abs(@mpwr-1)%8)>3
                        if (abs(@mpwr-1)%4)>1
                          zp4 = sqr(zp2)
                          z1p4 = sqr(z1p2)
                        else
                          zp4 = zri1^4
                          z1p4 = zjk1^4
                        endif
                      endif
                      if (abs(@mpwr-1)%16)>7
                        if (abs(@mpwr-1)%8)>3
                          zp8 = sqr(zp4)
                          z1p8 = sqr(z1p4)
                        elseif (abs(@mpwr-1)%4)>1
                          zp8 = zp2^4
                          z1p8 = z1p2^4
                        else
                          zp8 = sqr(zri1^4)
                          z1p8 = sqr(zjk1^4)
                        endif
                      endif
                      if abs(@mpwr-1)>15
                        if (abs(@mpwr-1)%16)>7
                          zp16 = sqr(zp8)
                          z1p16 = sqr(z1p8)
                        elseif (abs(@mpwr-1)%8)>3
                          zp16 = zp4^4
                          z1p16 = z1p4^4
                        elseif (abs(@mpwr-1)%4)>1
                          zp16 = sqr(zp2^4)
                          z1p16 = sqr(z1p2^4)
                        else
                          zp16 = (zri1^4)^4
                          z1p16 = (zjk1^4)^4
                        endif
                      endif
                      if abs(@mpwr-1)==31
                        zri1 = zp16*zp8*zp4*zp2*zri1
                        zjk1 = z1p16*z1p8*z1p4*z1p2*zjk1
                      elseif abs(@mpwr-1)==30
                        zri1 = zp16*zp8*zp4*zp2
                        zjk1 = z1p16*z1p8*z1p4*z1p2
                      elseif abs(@mpwr-1)==29
                        zri1 = zp16*zp8*zp4*zri1
                        zjk1 = z1p16*z1p8*z1p4*zjk1
                      elseif abs(@mpwr-1)==28
                        zri1 = zp16*zp8*zp4
                        zjk1 = z1p16*z1p8*z1p4
                      elseif abs(@mpwr-1)==27
                        zri1 = zp16*zp8*zp2*zri1
                        zjk1 = z1p16*z1p8*z1p2*zjk1
                      elseif abs(@mpwr-1)==26
                        zri1 = zp16*zp8*zp2
                        zjk1 = z1p16*z1p8*z1p2
                      elseif abs(@mpwr-1)==25
                        zri1 = zp16*zp8*zri1
                        zjk1 = z1p16*z1p8*zjk1
                      elseif abs(@mpwr-1)==24
                        zri1 = zp16*zp8
                        zjk1 = z1p16*z1p8
                      elseif abs(@mpwr-1)==23
                        zri1 = zp16*zp4*zp2*zri1
                        zjk1 = z1p16*z1p4*z1p2*zjk1
                      elseif abs(@mpwr-1)==22
                        zri1 = zp16*zp4*zp2
                        zjk1 = z1p16*z1p4*z1p2
                      elseif abs(@mpwr-1)==21
                        zri1 = zp16*zp4*zri1
                        zjk1 = z1p16*z1p4*zjk1
                      elseif abs(@mpwr-1)==20
                        zri1 = zp16*zp4
                        zjk1 = z1p16*z1p4
                      elseif abs(@mpwr-1)==19
                        zri1 = zp16*zp2*zri1
                        zjk1 = z1p16*z1p2*zjk1
                      elseif abs(@mpwr-1)==18
                        zri1 = zp16*zp2
                        zjk1 = z1p16*z1p2
                      elseif abs(@mpwr-1)==17
                        zri1 = zp16*zri1
                        zjk1 = z1p16*zjk1
                      elseif abs(@mpwr-1)==16
                        zri1 = zp16
                        zjk1 = z1p16
                      elseif abs(@mpwr-1)==15
                        zri1 = zp8*zp4*zp2*zri1
                        zjk1 = z1p8*z1p4*z1p2*zjk1
                      elseif abs(@mpwr-1)==14
                        zri1 = zp8*zp4*zp2
                        zjk1 = z1p8*z1p4*z1p2
                      elseif abs(@mpwr-1)==13
                        zri1 = zp8*zp4*zri1
                        zjk1 = z1p8*z1p4*zjk1
                      elseif abs(@mpwr-1)==12
                        zri1 = zp8*zp4
                        zjk1 = z1p8*z1p4
                      elseif abs(@mpwr-1)==11
                        zri1 = zp8*zp2*zri1
                        zjk1 = z1p8*z1p2*zjk1
                      elseif abs(@mpwr-1)==10
                        zri1 = zp8*zp2
                        zjk1 = z1p8*z1p2
                      elseif abs(@mpwr-1)==9
                        zri1 = zp8*zri1
                        zjk1 = z1p8*zjk1
                      elseif abs(@mpwr-1)==8
                        zri1 = zp8
                        zjk1 = z1p8
                      elseif abs(@mpwr-1)==7
                        zri1 = zp4*zp2*zri1
                        zjk1 = z1p4*z1p2*zjk1
                      elseif abs(@mpwr-1)==6
                        zri1 = zp4*zp2
                        zjk1 = z1p4*z1p2
                      elseif abs(@mpwr-1)==5
                        zri1 = zp4*zri1
                        zjk1 = z1p4*zjk1
                      elseif abs(@mpwr-1)==4
                        zri1 = zp4
                        zjk1 = z1p4
                      elseif abs(@mpwr-1)==3
                        zri1 = zp2*zri1
                        zjk1 = z1p2*zjk1
                      elseif abs(@mpwr-1)==2
                        zri1 = zp2
                        zjk1 = z1p2
                      endif
                    endif
                  endif
                else
                  zri1 = zri^(@mpwr-1.0)
                  zjk1 = zjk^(@mpwr-1.0)
                endif
                r = sqrt(th)
                zri = zri*zri1
                dzjk = r + flip(dzj)
                zjk = zjk*zjk1
                if r>0.0
                  dzri = zri1*dzri/r
                else
                  dzri = zri1
                endif
                if ph>0.0
                  dzjk = zjk1*dzjk/ph
                else
                  dzjk = zjk1
                endif
                magn = magn^@mpwr
                dzri = dr*real(dzjk)*dzri
                if @plus
                  dzj = dr*imag(dzjk)
                else
                  dzj = -dr*imag(dzjk)
                endif
                if @fractaltype==1
                  dzri = dzri + 1.0
                endif
                if @plus
                  zj = cj + magn*imag(zjk)
                else
                  zj = cj - magn*imag(zjk)
                endif
                zri = magn*real(zjk)*zri + cri
Logged

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

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
David Makin
Global Moderator
Fractal Senior
******
Posts: 2286



Makin' Magic Fractals
WWW
« Reply #32 on: January 09, 2010, 05:48:19 PM »

double th = atan2(y, x);
double ph = -atan2(z, sqrt(x * x + y * y));
x = r2 * cos(p * ph) * cos(p * th) + par.a;
y = r2 * cos(p * ph) * sin(p * th) + par.b;
z = r2 * sin(p * ph) + par.c;

Thanks - I realised it was just a matter of changing viewing angle and perspective to get the image so I removed the post asking which formula it was wink
Logged

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

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
David Makin
Global Moderator
Fractal Senior
******
Posts: 2286



Makin' Magic Fractals
WWW
« Reply #33 on: January 09, 2010, 06:23:55 PM »

How about timing? This 640*640 image took 35 seconds to render onscreen (so no anti-aliasing)
Using my program it takes 20 seconds to render 640x640 image (step distance was 0.3*DE).

Hi guys, just to point out that if you're comparing timings you really need to state the computing power used wink


Managed to get a similar view of the same Julia...


Rendered within UF @640*640 in 27secs on my 3GHz P4HT using the current optimised version of my formula, just to put that in context exactly the same image rendered using full trig for both the iterate and the derivative takes nearly 5 minutes.
Rendered using analytical DE, threshold 1e-3, step sizes = 0.45*DE.

Here's part of the optimised code:


Apologies - that code's only used for the Mandelbrot !
Here's part of the Julia version:

            elseif (r=|zri|)>0.0
              magn = sqrt(magn)
              if @usede
                dr = @mpwr*dr*magn^(@mpwr-1.0)
                if @fractaltype==1
                  dr = dr + 1.0
                endif
              endif
              if @mpwr==0.0
                zri = 1.0 + cri
                zj = cj
              elseif @mpwr==1.0
                zri = zri + cri
                if @plus
                  zj = cj + zj
                else
                  zj = cj - zj
                endif
              elseif @mpwr==2.0
                zri = sqr(zri)*(1.0 - sqr(zj)/r) + cri
                if @plus
                  zj = cj + 2.0*zj*sqrt(r)
                else
                  zj = cj - 2.0*zj*sqrt(r)
                endif
              elseif @mpwr==3.0
                zri = (zri^3)*(1.0 - 3.0*sqr(zj)/r) + cri
                if @plus
                  zj = cj - zj^3 + 3.0*zj*r
                else
                  zj = zj^3 - 3.0*zj*r + cj
                endif
              elseif @mpwr==4.0
                ztemp = ((r=sqrt(r)) + flip(zj))^4
                zri = real(ztemp)*(zri/r)^4 + cri
                if @plus
                  zj = cj + imag(ztemp)
                else
                  zj = cj - imag(ztemp)
                endif
              elseif @mpwr==round(@mpwr) && abs(@mpwr)<32
                zjk = (r=sqrt(r)) + flip(zj)
                if @mpwr<0.0
                  zri = r/zri
                  zjk = 1.0/zjk
                else
                  zri = zri/r
                endif
                if (abs(@mpwr)%4)>1
                  zp2 = sqr(zri)
                  z1p2 = sqr(zjk)
                endif
                if (abs(@mpwr)%8)>3
                  if (abs(@mpwr)%4)>1
                    zp4 = sqr(zp2)
                    z1p4 = sqr(z1p2)
                  else
                    zp4 = zri^4
                    z1p4 = zjk^4
                  endif
                endif
                if (abs(@mpwr)%16)>7
                  if (abs(@mpwr)%8)>3
                    zp8 = sqr(zp4)
                    z1p8 = sqr(z1p4)
                  elseif (abs(@mpwr)%4)>1
                    zp8 = zp2^4
                    z1p8 = z1p2^4
                  else
                    zp8 = sqr(zri^4)
                    z1p8 = sqr(zjk^4)
                  endif
                endif
                if abs(@mpwr)>15
                  if (abs(@mpwr)%16)>7
                    zp16 = sqr(zp8)
                    z1p16 = sqr(z1p8)
                  elseif (abs(@mpwr)%8)>3
                    zp16 = zp4^4
                    z1p16 = z1p4^4
                  elseif (abs(@mpwr)%4)>1
                    zp16 = sqr(zp2^4)
                    z1p16 = sqr(z1p2^4)
                  else
                    zp16 = (zri^4)^4
                    z1p16 = (zjk^4)^4
                  endif
                endif
                if abs(@mpwr)==31
                  zri = zp16*zp8*zp4*zp2*zri
                  zjk = z1p16*z1p8*z1p4*z1p2*zjk
                elseif abs(@mpwr)==30
                  zri = zp16*zp8*zp4*zp2
                  zjk = z1p16*z1p8*z1p4*z1p2
                elseif abs(@mpwr)==29
                  zri = zp16*zp8*zp4*zri
                  zjk = z1p16*z1p8*z1p4*zjk
                elseif abs(@mpwr)==28
                  zri = zp16*zp8*zp4
                  zjk = z1p16*z1p8*z1p4
                elseif abs(@mpwr)==27
                  zri = zp16*zp8*zp2*zri
                  zjk = z1p16*z1p8*z1p2*zjk
                elseif abs(@mpwr)==26
                  zri = zp16*zp8*zp2
                  zjk = z1p16*z1p8*z1p2
                elseif abs(@mpwr)==25
                  zri = zp16*zp8*zri
                  zjk = z1p16*z1p8*zjk
                elseif abs(@mpwr)==24
                  zri = zp16*zp8
                  zjk = z1p16*z1p8
                elseif abs(@mpwr)==23
                  zri = zp16*zp4*zp2*zri
                  zjk = z1p16*z1p4*z1p2*zjk
                elseif abs(@mpwr)==22
                  zri = zp16*zp4*zp2
                  zjk = z1p16*z1p4*z1p2
                elseif abs(@mpwr)==21
                  zri = zp16*zp4*zri
                  zjk = z1p16*z1p4*zjk
                elseif abs(@mpwr)==20
                  zri = zp16*zp4
                  zjk = z1p16*z1p4
                elseif abs(@mpwr)==19
                  zri = zp16*zp2*zri
                  zjk = z1p16*z1p2*zjk
                elseif abs(@mpwr)==18
                  zri = zp16*zp2
                  zjk = z1p16*z1p2
                elseif abs(@mpwr)==17
                  zri = zp16*zri
                  zjk = z1p16*zjk
                elseif abs(@mpwr)==16
                  zri = zp16
                  zjk = z1p16
                elseif abs(@mpwr)==15
                  zri = zp8*zp4*zp2*zri
                  zjk = z1p8*z1p4*z1p2*zjk
                elseif abs(@mpwr)==14
                  zri = zp8*zp4*zp2
                  zjk = z1p8*z1p4*z1p2
                elseif abs(@mpwr)==13
                  zri = zp8*zp4*zri
                  zjk = z1p8*z1p4*zjk
                elseif abs(@mpwr)==12
                  zri = zp8*zp4
                  zjk = z1p8*z1p4
                elseif abs(@mpwr)==11
                  zri = zp8*zp2*zri
                  zjk = z1p8*z1p2*zjk
                elseif abs(@mpwr)==10
                  zri = zp8*zp2
                  zjk = z1p8*z1p2
                elseif abs(@mpwr)==9
                  zri = zp8*zri
                  zjk = z1p8*zjk
                elseif abs(@mpwr)==8
                  zri = zp8
                  zjk = z1p8
                elseif abs(@mpwr)==7
                  zri = zp4*zp2*zri
                  zjk = z1p4*z1p2*zjk
                elseif abs(@mpwr)==6
                  zri = zp4*zp2
                  zjk = z1p4*z1p2
                elseif abs(@mpwr)==5
                  zri = zp4*zri
                  zjk = z1p4*zjk
                elseif abs(@mpwr)==4
                  zri = zp4
                  zjk = z1p4
                elseif abs(@mpwr)==3
                  zri = zp2*zri
                  zjk = z1p2*zjk
                elseif abs(@mpwr)==2
                  zri = zp2
                  zjk = z1p2
                endif
                zri = real(zjk)*zri + cri
                if @plus
                  zj = cj + imag(zjk)
                else
                  zj = cj - imag(zjk)
                endif
Logged

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

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
JosLeys
Strange Attractor
***
Posts: 258


WWW
« Reply #34 on: January 09, 2010, 07:37:39 PM »

Quote
If course  . My mistake. I don't know why I thought that R is the last value of triplex number after iterations. It should be:
drx=|R-Rx|/delta
dry=|R-Ry|/delta
drz=|R-Rz|/delta
It is exactly the same as you wrote. It means that your formulas are completely the same as mine.

Well, in the next step you are going to calculate sqrt(drx*drx+dry*dry+drz*drz), so the sign of R-Ri is of no importance..

Anyway, I tried your DE in UF, and it works, unfortunately whereas the image I posted earlier took 34 seconds, the one below took more than 12 minutes, so while it is a lot more accurate ( I see the spikes now also), it is a lot less efficient. I"m trying to find a way to optimise it..


* 4D_test_144a.jpg (432.4 KB, 634x634 - viewed 277 times.)
Logged
Buddhi
Fractal Iambus
***
Posts: 885



WWW
« Reply #35 on: January 10, 2010, 12:47:24 PM »

I have just finishing rendering z^2+c Mandelbulb in 8000x8000 resolution using delta DE algorithm. Step = 0.1*DE was used for very high quality. Image was downscaled to 4000x4000 for 2xAntialiasing


http://krzysztofmarczak.deviantart.com/art/Mandelbrot-3D-4000x4000-149897511
Logged

twinbee
Fractal Fertilizer
*****
Posts: 383



WWW
« Reply #36 on: January 10, 2010, 05:57:33 PM »

Wow, that's stunning Buddhi. I can even see some real fractal detail there too in certain places.
Logged
gaston3d
Guest
« Reply #37 on: January 10, 2010, 10:02:18 PM »

I have just finishing rendering z^2+c Mandelbulb in 8000x8000 resolution using delta DE algorithm. Step = 0.1*DE was used for very high quality. Image was downscaled to 4000x4000 for 2xAntialiasing

http://krzysztofmarczak.deviantart.com/art/Mandelbrot-3D-4000x4000-149897511

oak or maple leaf mandelbulb wink
Logged
fractalrebel
Fractal Lover
**
Posts: 211



WWW
« Reply #38 on: January 12, 2010, 01:37:25 AM »

Anyway, I tried your DE in UF, and it works, unfortunately whereas the image I posted earlier took 34 seconds, the one below took more than 12 minutes, so while it is a lot more accurate ( I see the spikes now also), it is a lot less efficient. I"m trying to find a way to optimise it..

I have sucessfully implemented Buddhi's DE in my UltraFractal formula. As Jos observed I get the little spike things sticking out. My default DE method (the Hart, Sandin and Kauffman method) gives a very similar (but not a good, imo wink) image. The render time was 1.58 min for a 640x640 image run on a DualCore E4500. I am still doing some optimization and testing, after which I will upload the modified formula to the UltraFractal database.
« Last Edit: January 12, 2010, 01:51:08 AM by fractalrebel » Logged

Pages: 1 2 [3]   Go Down
  Print  
 
Jump to:  

Related Topics
Subject Started by Replies Views Last post
Simple 2D IFS animation Movies Showcase (Rate My Movie) David Makin 0 806 Last post December 09, 2006, 11:41:48 PM
by David Makin
Plain and Simple Mandelbulb3D Gallery lenord 0 331 Last post February 27, 2011, 07:14:47 PM
by lenord
Non-trig expansion for RuckerBrot? Help & Support DarkBeam 1 831 Last post March 16, 2011, 01:41:42 PM
by DarkBeam
Does anyone know a trick to compute the derivative of Mandelbrot ? General Discussion tit_toinou 9 1789 Last post April 02, 2012, 03:53:12 AM
by David Makin
VTune Results; Compute vs. CVector3::IsNotANumber. Mandelbulber mancoast 0 435 Last post August 03, 2016, 03:42:16 AM
by mancoast

Powered by MySQL Powered by PHP Powered by SMF 1.1.21 | SMF © 2015, Simple Machines

Valid XHTML 1.0! Valid CSS! Dilber MC Theme by HarzeM
Page created in 0.351 seconds with 29 queries. (Pretty URLs adds 0.014s, 2q)