如何用C語言畫一個“聖誕樹”,我使用了左右鏡像的Sierpinski triangle,每層減去上方一小塊,再用符號點綴。可產生不同層數的「聖誕樹」,如下圖是5層的結果
#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的時候,就是只畫一主樹榦,樹榦越高就越幼:<img
n = 1的時候,利用遞迴畫向兩面分支,旋轉,越高的部分縮得越小。<img
n = 2 的時候,繼續分支出更細的樹支。n = 2 的時候,繼續分支出更細的樹支。<img
n = 3就差不多夠細節了。n = 3就差不多夠細節了。
代碼長一點,為了容易理解我不「壓縮」它了。
#include <math.h>#include <stdio.h>#include <stdlib.h>#define PI 3.14159265359float 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: 簡單地加入裝飾及絲帶,在命令列可以選擇放大倍率,下圖是兩倍大的。
<img src="https://pic2.zhimg.com/fa09e223f37b214d5bca14953366150d_b.jpg" data-rawwidth="711" data-rawheight="823" class="origin_image zh-lightbox-thumb" width="711" data-original="https://pic2.zhimg.com/fa09e223f37b214d5bca14953366150d_r.jpg">// 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.14159265359ffloat 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)))); } } // ...}
其實我在回答這問題的時候,並沒有計劃,只是一步一步地嘗試。現在我覺得用這規模的代碼大概不能再怎麼進展了。不過今天看到大堂裡的聖誕樹,覺得那些裝飾物還挻有趣的,有時候除了畫整體,也可以畫局部,看看是否能再更新。
聖誕節快樂!