渲染mandelbrot集的最快算法是什么? [英] What are the fastest algorithms for rendering the mandelbrot set?

查看:57
本文介绍了渲染mandelbrot集的最快算法是什么?的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我尝试了许多用于渲染Mandelbrot集的算法,其中包括朴素的逸出时间算法以及优化的逸出时间算法.但是,是否有更快的算法可以像我们在YouTube上看到的那样有效地产生真正的深度缩放.另外,我很想了解如何提高C/C ++ double

之外的精度

解决方案

即使是高端CPU也要比普通GPU慢得多.即使使用GPU上的幼稚迭代算法,您也可以进行实时渲染.因此,在GPU上使用更好的算法可以实现较高的缩放比例,但是对于您需要的任何体面算法,都是如此:

  • 多次渲染,因为我们无法在GPU上自行修改纹理
  • 高精度浮点,因为 floats/doubles 是不够的.

相关的质量检查很少:

  • 因此它以64x32 Mandelbrot像素显示为128x64抖动图像.在我的AVR32上,这可能比幼稚的方法(可能是3-4fps)快8倍...代码可能会得到更优化,但是请记住,Mandelbrot并不是唯一运行的东西,因为我在后台有一些ISR处理程序来处理LCD以及我的

    RT GIF CPU捕获(n = 62 ++,GIF 4倍缩小):

    I've tried many algorithms for the rendering of the Mandelbrot set, inclusive of the naive escape time algorithm, as well as the optimized escape time algorithm. But, are there faster algorithms that are used to produce really deep zooms efficiently like the ones we see on YouTube. Also, I would love to get some ideas on how to increase my precision beyond the C/C++ double

    解决方案

    Even High end CPU will be much slower in comparison to average GPU. You can get to real time rendering even with naive iteration algo on GPU. So using better algorithms on GPU could get to high zooms however for any decent algo you need:

    • multi pass rendering as we can not self modify a texture on GPU
    • high precision floating point as floats/doubles are not enough.

    Here few related QAs:

    which might get you kick started...

    One way to speed up is you can use fractional escape like I did in the first link. It improves image quality while keeping max iteration low.

    The second link will get you approximation of which parts of fractal are in and out and how far. Its not very accurate but can be used to avoid computing iterations for parts that are "outside for sure".

    Next link will show you how to achieve better precision.

    Last link is about Perturbation The idea is that you use high precision math only for some reference points and use that to compute its neighbor points with low precision math without loosing precision. Never used that however but looks promising.

    And finally once you achieved fast rendering you might want to aim for this:

    Here a small example of 3* 64bit double used for single value in GLSL:

    // high precision float (very slow)
    dvec3 fnor(dvec3 a)
        {
        dvec3 c=a;
        if (abs(c.x)>1e-5){ c.y+=c.x; c.x=0.0; }
        if (abs(c.y)>1e+5){ c.z+=c.y; c.y=0.0; }
        return c;
        }
    double fget(dvec3 a){ return a.x+a.y+a.z; }
    dvec3 fset(double a){ return fnor(dvec3(a,0.0,0.0)); }
    dvec3 fadd(dvec3 a,double b){ return fnor(a+fset(b)); }
    dvec3 fsub(dvec3 a,double b){ return fnor(a-fset(b)); }
    dvec3 fmul(dvec3 a,double b){ return fnor(a*b); }
    dvec3 fadd(dvec3 a,dvec3 b){ return fnor(a+b); }
    dvec3 fsub(dvec3 a,dvec3 b){ return fnor(a-b); }
    dvec3 fmul(dvec3 a,dvec3 b)
        {
        dvec3 c;
        c =fnor(a*b.x);
        c+=fnor(a*b.y);
        c+=fnor(a*b.z);
        return fnor(c);
        }
    

    so each hi precision value is dvec3 ... the thresholds in fnor can be changed to any ranges. You can convert this to vec3 and float ...

    [Edit1] "fast" C++ example

    Ok I wanted to try my new SSD1306 driver along with my AVR32 MCU to compute Mandelbrot So I can compare speed with this Arduino + 3D + Pong + Mandelbrot. I used AT32UC3A3256 with ~66MHz no FPU no GPU and 128x64x1bpp display. No external memory only internal 16+32+32 KByte. Naive Mandlebrot was way to slow (~2.5sec per frame) so I busted up something like this (taking advantage of that position and zoom of the view is sort of continuous):

    1. reduce resolution by 2

      to make room for dithering as my output is just B&W

    2. use variable max iteration n based on zoom

      On change of n invalidate last frame to enforce full recompute. I know this is slow but it happens only 3 times on transitions between zoom ranges.

      Scaling count from last frame is not looking good as its not linear.

      Its possible to use the last counts but for that it would be needed also to remember the complex variables used for iteration and that would take too much memory.

    3. remember last frame and also which x,y screen coordinate mapped to which Mandelbrot coordinate.

    4. On each frame compute the mapping between screen coordinates and Mandelbrot coordinates.

    5. remap last frame to adjust to new position and zoom

      so simply look at the data from #3,#4 and if we have there the same positions in both last and actual frame (closer then half of pixel size), copy the pixels. And recompute the rest.

      This will hugely improve performance if your view is smooth (so position and zoom does not change a lot on per frame basis).

    I know its a bit vague description so here a C++ code where you can infer all doubts:

    //---------------------------------------------------------------------------
    //--- Fast Mandelbrot set ver: 1.000 ----------------------------------------
    //---------------------------------------------------------------------------
    template<int xs,int ys,int sh> void mandelbrot_draw(float mx,float my,float zoom)
        {
        // xs,ys - screen resolution
        // sh    - log2(pixel_size) ... dithering pixel size
        // mx,my - Mandelbrot position (center of view) <-1.5,+0.5>,<-1.0,+1.0>
        // zoom  - zoom
        // ----------------
        // (previous/actual) frame
        static U8 p[xs>>sh][ys>>sh];            // intensities (raw Mandelbrot image)
        static int n0=0;                        // max iteraions
        static float px[(xs>>sh)+1]={-1000.0};  // pixel x position in Mandlebrot
        static float py[(ys>>sh)+1];            // pixel y position in Mandlebrot
        // temp variables
        U8 shd;                                 // just pattern for dithering
        int ix,iy,i,n,jx,jy,kx,ky,sz;           // index variables
        int nx=xs>>sh,ny=ys>>sh;                // real Mandelbrot resolution
        float fx,fy,fd;                         // floating Mandlebrot position and pixel step
        float x,y,xx,yy,q;                      // Mandelbrot iteration stuff (this need to be high precision)
        int   qx[xs>>sh],qy[ys>>sh];            // maping of pixels between last and actual frame
        float px0[xs>>sh],py0[ys>>sh];          // pixel position in Mandlebrot from last frame
        // init vars
             if (zoom<   10.0) n= 31;
        else if (zoom<  100.0) n= 63;
        else if (zoom< 1000.0) n=127;
        else                   n=255;
        sz=1<<sh;
        ix=xs; if (ix>ys) ix=ys; ix/=sz;
        fd=2.0/(float(ix-1)*zoom);
        mx-=float(xs>>(1+sh))*fd;
        my-=float(ys>>(1+sh))*fd;
        // init buffers
        if ((px[0]<-999.0)||(n0!=n))
            {
            n0=n;
            for (ix=0;ix<nx;ix++) px[ix]=-999.0; 
            for (iy=0;iy<ny;iy++) py[iy]=-999.0;
            for (ix=0;ix<nx;ix++)
             for (iy=0;iy<ny;iy++)
              p[ix][iy]=0;
            }
        // store old and compute new float positions of pixels in Mandelbrot to px[],py[],px0[],py0[]
        for (fx=mx,ix=0;ix<nx;ix++,fx+=fd){ px0[ix]=px[ix]; px[ix]=fx; qx[ix]=-1; }
        for (fy=my,iy=0;iy<ny;iy++,fy+=fd){ py0[iy]=py[iy]; py[iy]=fy; qy[iy]=-1; }
        // match old and new x coordinates to qx[]
        for (ix=0,jx=0;(ix<nx)&&(jx<nx);)
            {
            x=px[ix]; y=px0[jx];
            xx=(x-y)/fd; if (xx<0.0) xx=-xx;
            if (xx<=0.5){ qx[ix]=jx; px[ix]=y; }
            if (x<y) ix++; else jx++;
            }
        // match old and new y coordinates to qy[]
        for (ix=0,jx=0;(ix<ny)&&(jx<ny);)
            {
            x=py[ix]; y=py0[jx];
            xx=(x-y)/fd; if (xx<0.0) xx=-xx;
            if (xx<=0.5){ qy[ix]=jx; py[ix]=y; }
            if (x<y) ix++; else jx++;
            }
        // remap p[][] by qx[]
        for (ix=0,jx=nx-1;ix<nx;ix++,jx--)
            {
            i=qx[ix]; if ((i>=0)&&(i>=ix)) for (iy=0;iy<ny;iy++) p[ix][iy]=p[i][iy];
            i=qx[jx]; if ((i>=0)&&(i<=jx)) for (iy=0;iy<ny;iy++) p[jx][iy]=p[i][iy];
            }
        // remap p[][] by qy[]
        for (iy=0,jy=ny-1;iy<ny;iy++,jy--)
            {
            i=qy[iy]; if ((i>=0)&&(i>=iy)) for (ix=0;ix<nx;ix++) p[ix][iy]=p[ix][i];
            i=qy[jy]; if ((i>=0)&&(i<=jy)) for (ix=0;ix<nx;ix++) p[ix][jy]=p[ix][i];
            }
        // Mandelbrot
        for (iy=0,ky=0,fy=py[iy];iy<ny;iy++,ky+=sz,fy=py[iy]) if ((fy>=-1.0)&&(fy<=+1.0))
         for (ix=0,kx=0,fx=px[ix];ix<nx;ix++,kx+=sz,fx=px[ix]) if ((fx>=-1.5)&&(fx<=+0.5))
            {
            // invalid qx,qy ... recompute Mandelbrot
            if ((qx[ix]<0)||(qy[iy]<0))
                {
                for (x=0.0,y=0.0,xx=0.0,yy=0.0,i=0;(i<n)&&(xx+yy<4.0);i++)
                    {
                    q=xx-yy+fx;
                    y=(2.0*x*y)+fy;
                    x=q;
                    xx=x*x;
                    yy=y*y;
                    }
                i=(16*i)/(n-1); if (i>16) i=16; if (i<0) i=0;
                i=16-i; p[ix][iy]=i;
                }
            // use stored intensity
            else i=p[ix][iy];
            // render point with intensity i coresponding to ix,iy position in map
            for (i<<=3                 ,jy=0;jy<sz;jy++)
             for (shd=shade8x8[i+(jy&7)],jx=0;jx<sz;jx++)
              lcd.pixel(kx+jx,ky+jy,shd&(1<<(jx&7)));
            }
        }
    //---------------------------------------------------------------------------
    //---------------------------------------------------------------------------
    //---------------------------------------------------------------------------
    

    the lcd and shade8x8 stuff can be found in the linked SSD1306 QA. However you can ignore it its just dithering and outputting a pixel so you can instead output the i directly (even without the scaling to <0..16>.

    Here preview (on PC as I was lazy to connect camera ...):

    so it 64x32 Mandelbrot pixels displayed as 128x64 dithered image. On my AVR32 is this may be 8x time faster than naive method (maybe 3-4fps)... The code might be more optimized however take in mind Mandelbrot is not the only stuff running as I have some ISR handlers on backround to handle the LCD and also my TTS engine based on this which I upgraded a lot since then and use for debugging of this (yes it can speek in parallel to rendering). Also I am low on memory as my 3D engine takes a lot ~11 KByte away (mostly depth buffer).

    The preview was done with this code (inside timer):

    static float zoom=1.0;
    mandelbrot_draw<128,64,1>(+0.37,-0.1,zoom);
    zoom*=1.02; if (zoom>100000) zoom=1.0;
    

    Also for non AVR32 C++ environment use this:

    //------------------------------------------------------------------------------------------
    #ifndef _AVR32_compiler_h
    #define _AVR32_compiler_h
    //------------------------------------------------------------------------------------------
    typedef int32_t  S32;
    typedef int16_t  S16;
    typedef int8_t   S8;
    typedef uint32_t U32;
    typedef uint16_t U16;
    typedef uint8_t  U8;
    //------------------------------------------------------------------------------------------
    #endif
    //------------------------------------------------------------------------------------------
    

    [Edit2] higher float precision in GLSL

    The main problem with Mandelbrot is that it need to add numbers with very big exponent difference. For +,- operations we need to align mantissa of both operands, add them as integer and normalize back to scientific notation. However if the exponent difference is big then the result mantissa need more bits than can fit into 32 bit float so only 24 most significant bits are preserved. This creates the rounding errors causing your pixelation. If you look at 32bit float in binary you will see this:

    float a=1000.000,b=0.00001,c=a+b; 
    
    //012345678901234567890123456789 ... just to easy count bits
    a=1111101000b                                         // a=1000
    b=         0.00000000000000001010011111000101101011b  // b=0.00000999999974737875
    c=1111101000.00000000000000001010011111000101101011b  // not rounded result
    c=1111101000.00000000000000b                          // c=1000 rounded to 24 bits of mantissa
    

    Now the idea is to enlarge the number of mantissa bits. The easiest trick is to have 2 floats instead of one:

    //012345678901234567890123 ... just to easy count bits
    a=1111101000b                                         // a=1000
                           //012345678901234567890123 ... just to easy count     b=         0.0000000000000000101001111100010110101100b    // b=0.00000999999974737875
    c=1111101000.00000000000000001010011111000101101011b  // not rounded result
    c=1111101000.00000000000000b                          // c=1000 rounded to 24 
     +          .0000000000000000101001111100010110101100b
                               //012345678901234567890123 ... just to easy count bits
    

    so some part of the result is in one float and the rest in the other... The more floats per single value we have the bigger the mantissa. However doing this on bit exact division of big mantissa into 24 bit chunks would be complicated and slow in GLSL (if even possible due to GLSL limitations). Instead we can select for each of the float some range of exponents (just like in example above).

    So in the example we got 3 floats (vec3) per single (float) value. Each of the coordinates represent different range:

           abs(x) <= 1e-5
    1e-5 < abs(y) <= 1e+5
    1e+5 < abs(z) 
    

    and value = (x+y+z) so we kind of have 3*24 bit mantissa however the ranges are not exactly matching 24 bits. For that the exponent range should be divided by:

    log10(2^24)=7.2247198959355486851297334733878
    

    instead of 10... for example something like this:

           abs(x) <= 1e-7
    1e-7 < abs(y) <= 1e+0
    1e+0 < abs(z) 
    

    Also the ranges must be selected so they handle the ranges of values you use otherwise it would be for nothing. So if your numbers are <4 its pointless to have range >10^+5 So first you need to see what bounds of values you have, then dissect it to exponent ranges (as many as you have floats per value).

    Beware some (but much less than native float) rounding still occurs !!!

    Now doing operations on such numbers is slightly more complicated than on normal floats as you need to handle each value as bracketed sum of all components so:

    (x0+y0+z0) + (x1+y1+z1) = (x0+x1 + y0+y1 + z0+z1)
    (x0+y0+z0) - (x1+y1+z1) = (x0-x1 + y0-y1 + z0-z1)
    (x0+y0+z0) * (x1+y1+z1) = x0*(x1+y1+z1) + y0*(x1+y1+z1) + z0*(x1+y1+z1)
    

    And do not forget to normalize the values back to the defined ranges. Avoid adding small and big (abs) values so avoid x0+z0 etc ...

    [Edit3] new win32 demo CPU vs. GPU

    Both executables are preset to the same location and zoom to show when the doubles starts to round off. I had to upgrade slightly the way how px,py coordinates are computed as around 10^9 the y axis started to deviate at this location (The threshold still might be too big for other locations)

    Here preview CPU vs. GPU for high zoom (n=1600):

    RT GIF capture of CPU (n=62++, GIF 4x scaled down):

    这篇关于渲染mandelbrot集的最快算法是什么?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆