2009年10月31日 星期六

Fractal

 
底下秀出 4 種圖:

1.Koch curve
2.Sierpinski’s Gasket
3.Tree curve
4.Julia set

其中,sierpinski 外緣線需自行接上,
Julia set 觀察窗需手動調整:

#include <windows.h>
#include <math.h>
#include <complex>

#define RGB(r,g,b) (b<<16)|(g<<8)|r
enum {W=800, H=600};
HDC hdc;
HANDLE  hwnd, hIn, hOut;

namespace Turtle
{
    //--------------------- 作圖 -----------------------

    const float rd = 3.1415927/180.0;
    int* pic, X, Y, A, ang, step, C;   //定義: 畫布, 座標, 角度, 步距, 顏色
    float scale;
    float branch;

    void clear ()
    {
        int i = W*H; while (i--) pic[i] = 0;
    }
    void show ()                //注意: 顯示器色彩格式需為 32-bits
    {
        BITMAPINFO info = {{40, W, -H, 1,32,0,0,0,0,0,0}, {{0}}};
        SetDIBitsToDevice 
           (hdc, 0,0,W,H, 0,0,0,H, pic, &info, DIB_RGB_COLORS);
    }
    void set (int x, int y, int d=0, int a=0) {X=x; Y=y; ang=a; step=d;}
    void setxy (int x, int y) {set(x,y);}
    void pixel (int x, int y, int c)                //繪點
    {
        if (x<0 || y<0 || x>=W || y>=H) return;
        pic[y*W+x] = c;
    }  
    void line (int x, int y, int x2, int y2, int c=0xffffff)
    {
        int dx= x2-x, dy= y2-y, dirx=1, diry=1, b=0;
        if (dx < 0)  {dx = -dx; dirx= -1;}           
        if (dy < 0)  {dy = -dy; diry= -1;}           
        pixel (x,y,c);
        if (dx<=1 && dy<=1)  return;                 //兩點間隔太小,不予處理
        if (dx < dy) {                               //以Y軸為基準繪製線段
            int hy= dy>>1, dotn= dy-2, n=1;
            do {b += dx;   if (b > hy)  {b -= dy; x += dirx;}
                y += diry; pixel (x,y,c); n++;
            }while (dotn--);
        }else{                                       //以X軸為基準繪製線段
            int hx= dx>>1, dotn= dx-2, n=1;
            do {b += dy;   if (b > hx)  {b -= dx; y += diry;}
                x += dirx; pixel (x,y,c); n++;
            }while (dotn--);
        }
    }
    void move (int n, int c=0xffffff)
    {
        int x = n* cos (rd*A), y = n* sin (rd*A);
        line (X, Y, X+x, Y+y, c);
        X+=x; Y+=y;
    }
    void mov (int c=0xffffff) {move (step,c);}
    void turn (int a)
    {
        A+=a;
        if (A>360) A-=360;
        else if (A<360) A+=360;
    }

    //--------------------- 碎形 -----------------------

    void koch (int n)
    {
        if (0==n) {mov(); return;}
        koch (n-1);  turn (-60);
        koch (n-1);  turn (120);
        koch (n-1);  turn (-60);
        koch (n-1);
    }
    void sierpinski (int n)
    {
        if (0==n) {mov(); return;}
        sierpinski (n-1); turn(-120);
        sierpinski (n-1); turn(120);
        sierpinski (n-1); turn(120);
        sierpinski (n-1); turn(-120);
        sierpinski (n-1);
    }
    void tree (int n, int x, int y, float len, float ang)
    {
        if (0==n) return;
        setxy(x,y); A=-ang;  move(len); x=X; y=Y;
        tree (n-1, x,y, len/scale, ang-branch);
        tree (n-1, x,y, len/scale, ang+branch);
    }
    void julia (int n)
    {
        // z = x+iy,  F(z) = z^2+C,  |F(z)| = x^2+y^2
        int i, k;
        float v,xp,yp, x0,y0,x1,y1,x,y,z,cx,cy;
        cx=-0.5;  cy=0.5;
        for (x0=-2.0; x0<=2.0; x0+=0.005)
            for (y0=-2.0; y0<=2.0; y0+=0.005){
                for (x = x0, y = y0, k = 0;
                    x*x+y*y<4.0 && k<n; ++k){    //迭代 n 次,直到數列發散
                    x1 = x*x - y*y + cx;
                    y1 = 2*x*y + cy;
                    x = x1; y = y1;
                }
                xp = W/4*x0+W/2;
                yp = -H/4*y0+H/2;
                v = x*x + y*y;
                if (v < 4.0) {
                    i = v*64;
                    pixel (xp, yp, RGB (100,220,i));
                }else
                    pixel (xp, yp, RGB (30,100,200));
            }
    }
};
using namespace Turtle;


int main (/* daviddr 081020 */)
{
    SetConsoleTitle (L"Conso");
    hdc  = GetDC (FindWindowW (0, L"Conso"));
    hIn  = GetStdHandle (STD_INPUT_HANDLE);
    hOut = GetStdHandle (STD_OUTPUT_HANDLE);
    SMALL_RECT size = {0, 0, W-1, H-1};      
    SetConsoleWindowInfo (hOut, TRUE, &size);   //變更視窗大小
    Turtle::pic = new int[W*H];    //畫布
   
    clear ();
    set (20,380, 8);  koch (4);   
    show ();
    system ("pause");

    clear ();
    set (20,380, 8);  sierpinski (5);   
    show ();
    system ("pause");

    clear ();
    scale = 1.3;
    branch = 23;
    tree (8, 200,380, 64, 80);   
    show ();
    system ("pause");

    clear();
    julia(32);
    show ();
    system ("pause");

    delete [] Turtle::pic;
    DeleteDC (hdc);
    CloseHandle (hIn);
    CloseHandle (hOut);

    return 0;
}
sierpinski 外緣線由小而大逐次爬行的作法,
很難排除累計誤差,無法形成正三角。
可改為由外而內,用內插方式來算落點:

private:
    void sierpinski (int n, int x1, int y1, int x2, int y2, int x3, int y3)   
    {
        if (0 == n) {            
            setxy (x1, y1);
            mov(); turn (-120);
            mov(); turn (-120);
            mov(); turn (-120);
            return;
        }
        sierpinski (n-1, x1,y1, x3,y1, (x1+x3)/2, (y1+y3)/2);
        sierpinski (n-1, x3,y1, x2,y2, (x2+x3)/2, (y2+y3)/2);
        sierpinski (n-1, (x1+x3)/2, (y1+y3)/2, (x2+x3)/2, (y2+y3)/2, x3, y3);
    }

public:
    void sierpin2 (int n)                       //自頂至底
    {
        int d = step * (1<<n);
        sierpinski (n, X,Y, X+d,Y, X+d/2, int(float(Y)-1.72384*d/2));
    }

下面是呼叫 koch (4), sierpin2 (5), tree (8,..) 的結果:


沒有留言:

張貼留言