如何用 C 語言畫一個“聖誕樹”？

知乎 - Milo Yip發表於2015-12-25

我使用了左右映象的Sierpinski triangle，每層減去上方一小塊，再用符號點綴。可生成不同層數的「聖誕樹」，如下圖是5層的結果。

#include <stdio.h>
#include <stdlib.h>

int main(int argc, char* argv[]) {
    int n = argc > 1 ? atoi(argv[1]) : 4;
    for (int j = 1; j <= n; j++) {
        int s = 1 << j, k = (1 << n) - s, x;
        for (int y = s - j; y >= 0; y--, putchar('\n')) {
            for (x = 0; x < y + k; x++) printf("  ");
            for (x = 0; x + y < s; x++) printf("%c ", '!' ^ y & x);
            for (x = 1; x + y < s; x++) printf("%c ", '!' ^ y & (s - y - x - 1));
        }
    }
}

基本程式碼來自Sierpinski triangle的實現，字元的想法來自於code golf – Draw A Sierpinski Triangle。
–
更新1: 上面的是我嘗試儘量用最少程式碼來畫一個抽象一點的聖誕樹，因此樹幹都沒有。然後，我嘗試用更真實一點的風格。因為樹是一個比較自相似的形狀，這次使用遞迴方式描述樹幹和分支。

n = 0的時候，就是隻畫一主樹幹，樹幹越高就越幼：

n = 1的時候，利用遞迴畫向兩面分支，旋轉，越高的部分縮得越小。

n = 2 的時候，繼續分支出更細的樹支。

n = 3就差不多夠細節了。

程式碼長一點，為了容易理解我不「壓縮」它了。

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define PI 3.14159265359

float sx, sy;

float sdCircle(float px, float py, float r) {
    float dx = px - sx, dy = py - sy;
    return sqrtf(dx * dx + dy * dy) - r;
}

float opUnion(float d1, float d2) {
    return d1 < d2 ? d1 : d2;
}

#define T px + scale * r * cosf(theta), py + scale * r * sin(theta)

float f(float px, float py, float theta, float scale, int n) {
    float d = 0.0f;
    for (float r = 0.0f; r < 0.8f; r += 0.02f)
        d = opUnion(d, sdCircle(T, 0.05f * scale * (0.95f - r)));

    if (n > 0)
        for (int t = -1; t <= 1; t += 2) {
            float tt = theta + t * 1.8f;
            float ss = scale * 0.9f;
            for (float r = 0.2f; r < 0.8f; r += 0.1f) {
                d = opUnion(d, f(T, tt, ss * 0.5f, n - 1));
                ss *= 0.8f;
            }
        }

    return d;
}

int main(int argc, char* argv[]) {
    int n = argc > 1 ? atoi(argv[1]) : 3;
    for (sy = 0.8f; sy > 0.0f; sy -= 0.02f, putchar('\n'))
        for (sx = -0.35f; sx < 0.35f; sx += 0.01f)
            putchar(f(0, 0, PI * 0.5f, 1.0f, n) < 0 ? '*' : ' ');
}

這段程式碼實際上是用了圓形的距離場來建模，並且沒有優化。

這是一棵「祼樹」，未能稱得上是「聖誕樹」。

–
更新2: 簡單地加入裝飾及絲帶，在命令列可以選擇放大倍率，下圖是兩倍大的。

// f() 及之前的部分沿上

int ribbon() {
    float x = (fmodf(sy, 0.1f) / 0.1f - 0.5f) * 0.5f;
    return sx >= x - 0.05f && sx <= x + 0.05f;
}

int main(int argc, char* argv[]) {
    int n = argc > 1 ? atoi(argv[1]) : 3;
    float zoom = argc > 2 ? atof(argv[2]) : 1.0f;
    for (sy = 0.8f; sy > 0.0f; sy -= 0.02f / zoom, putchar('\n'))
        for (sx = -0.35f; sx < 0.35f; sx += 0.01f / zoom) {
            if (f(0, 0, PI * 0.5f, 1.0f, n) < 0.0f) {
                if (sy < 0.1f)
                    putchar('.');
                else {
                    if (ribbon())
                        putchar('=');
                    else
                        putchar("............................#j&o"[rand() % 32]);
                }
            }
            else
                putchar(' ');
        }
}

2D的我想已差不多了。接下來看看有沒有空嘗試3D的。

–
更新3：終於要3D了。之前每個節點是往左和右分支，在三維中我們可以更自由一點，我嘗試在每個節點申出6個分支。最後用了簡單的Lambertian著色（即max(dot(N, L), 0）。

n = 1 的時候比較容易看出立體的著色：

可是n=3的時候已亂得難以辨認：

估計是因為aliasing而做成的。由於光照已經使用了finite difference來計演算法線，效能已經很差，我就不再嘗試做Supersampling去解決aliasing的問題了。另外也許Ambient occlusion對這問題也有幫助，不過需要更多的取樣。

因為需要三維旋轉，不能像二維簡單使用一個角度來代表旋轉，所以這段程式碼加入了不少矩陣運算。當然用四元數也是可以的。

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define PI 3.14159265359f

float sx, sy;

typedef float Mat[4][4];
typedef float Vec[4];

void scale(Mat* m, float s) {
    Mat temp = { {s,0,0,0}, {0,s,0,0 }, { 0,0,s,0 }, { 0,0,0,1 } };
    memcpy(m, &temp, sizeof(Mat));
}

void rotateY(Mat* m, float t) {
    float c = cosf(t), s = sinf(t);
    Mat temp = { {c,0,s,0}, {0,1,0,0}, {-s,0,c,0}, {0,0,0,1} };
    memcpy(m, &temp, sizeof(Mat));
}

void rotateZ(Mat* m, float t) {
    float c = cosf(t), s = sinf(t);
    Mat temp = { {c,-s,0,0}, {s,c,0,0}, {0,0,1,0}, {0,0,0,1} };
    memcpy(m, &temp, sizeof(Mat));
}

void translate(Mat* m, float x, float y, float z) {
    Mat temp = { {1,0,0,x}, {0,1,0,y}, {0,0,1,z}, {0,0,0,1} };
    memcpy(m, &temp, sizeof(Mat));
}

void mul(Mat* m, Mat a, Mat b) {
    Mat temp;
    for (int j = 0; j < 4; j++)
        for (int i = 0; i < 4; i++) {
            temp[j][i] = 0.0f;
            for (int k = 0; k < 4; k++)
                temp[j][i] += a[j][k] * b[k][i];
        }
    memcpy(m, &temp, sizeof(Mat));    
}

void transformPosition(Vec* r, Mat m, Vec v) {
    Vec temp = { 0, 0, 0, 0 };
    for (int j = 0; j < 4; j++)
        for (int i = 0; i < 4; i++)
            temp[j] += m[j][i] * v[i];
    memcpy(r, &temp, sizeof(Vec));    
}

float transformLength(Mat m, float r) {
    return sqrtf(m[0][0] * m[0][0] + m[0][1] * m[0][1] + m[0][2] * m[0][2]) * r;
}

float sphere(Vec c, float r) {
    float dx = c[0] - sx, dy = c[1] - sy;
    float a = dx * dx + dy * dy;
    return a < r * r ? sqrtf(r * r - a) + c[2] : -1.0f;
}

float opUnion(float z1, float z2) {
    return z1 > z2 ? z1 : z2;
}

float f(Mat m, int n) {
    float z = -1.0f;
    for (float r = 0.0f; r < 0.8f; r += 0.02f) {
        Vec v = { 0.0f, r, 0.0f, 1.0f };
        transformPosition(&v, m, v);
        z = opUnion(z, sphere(v, transformLength(m, 0.05f * (0.95f - r))));
    }

    if (n > 0) {
        Mat ry, rz, s, t, m2, m3;
        rotateZ(&rz, 1.8f);

        for (int p = 0; p < 6; p++) {
            rotateY(&ry, p * (2 * PI / 6));
            mul(&m2, ry, rz);
            float ss = 0.45f;
            for (float r = 0.2f; r < 0.8f; r += 0.1f) {
                scale(&s, ss);
                translate(&t, 0.0f, r, 0.0f);
                mul(&m3, s, m2);
                mul(&m3, t, m3);
                mul(&m3, m, m3);
                z = opUnion(z, f(m3, n - 1));
                ss *= 0.8f;
            }
        }
    }

    return z;
}

float f0(float x, float y, int n) {
    sx = x;
    sy = y;
    Mat m;
    scale(&m, 1.0f);
    return f(m, n);
}

int main(int argc, char* argv[]) {
    int n = argc > 1 ? atoi(argv[1]) : 3;
    float zoom = argc > 2 ? atof(argv[2]) : 1.0f;
    for (float y = 0.8f; y > -0.0f; y -= 0.02f / zoom, putchar('\n'))
        for (float x = -0.35f; x < 0.35f; x += 0.01f / zoom) {
            float z = f0(x, y, n);
            if (z > -1.0f) {
                float nz = 0.001f;
                float nx = f0(x + nz, y, n) - z;
                float ny = f0(x, y + nz, n) - z;
                float nd = sqrtf(nx * nx + ny * ny + nz * nz);
                float d = (nx - ny + nz) / sqrtf(3) / nd;
                d = d > 0.0f ? d : 0.0f;
                // d = d < 1.0f ? d : 1.0f;
                putchar(".-:=+*#%@@"[(int)(d * 9.0f)]);
            }
            else
                putchar(' ');
        }
}

更新4：發現之前的TransformLength()寫錯了，上面已更正。另外，考慮提升效能時，一般是需要一些空間剖分的方式去加速檢查，但這裡剛好是一個樹狀的場景結構，可以簡單使用Bounding volume hierarchy，我使用了球體作為包圍體積。只需加幾句程式碼，便可以大大縮減執行時間。

另外，考慮到太小的葉片是很難取樣得到好看的結果，我嘗試以一個較大的球體去表現葉片（就如素描時考慮更整體的光暗而不是每片葉片的光暗），我覺得結果有進步。

float f(Mat m, int n) {
    // Culling
    {
        Vec v = { 0.0f, 0.5f, 0.0f, 1.0f };
        transformPosition(&v, m, v);        
        if (sphere(v, transformLength(m, 0.55f)) == -1.0f)
            return -1.0f;
    }

    float z = -1.0f;

    if (n == 0) { // Leaf
        Vec v = { 0.0f, 0.5f, 0.0f, 1.0f };
        transformPosition(&v, m, v);        
        z = sphere(v, transformLength(m, 0.3f));
    }  
    else { // Branch
        for (float r = 0.0f; r < 0.8f; r += 0.02f) {
            Vec v = { 0.0f, r, 0.0f, 1.0f };
            transformPosition(&v, m, v);
            z = opUnion(z, sphere(v, transformLength(m, 0.05f * (0.95f - r))));
        }
    }

    // ...
}

其實我在回答這問題的時候，並沒有計劃，只是一步一步地嘗試。現在我覺得用這規模的程式碼大概不能再怎麼進展了。不過今天看到大堂裡的聖誕樹，覺得那些裝飾物還挻有趣的，有時候除了畫整體，也可以畫區域性，看看是否能再更新。

相關文章