用C语言画希尔伯特曲线

1,342 阅读5分钟
原文链接: zhuanlan.zhihu.com

希尔伯特曲线Hilbert curve)是一种空间填充曲线(space filling curve)的分形图案,由德国数学家 David Hilbert(1862-1943)发现 [1]。简单来说,空间填充曲线就是可以用「一笔画」的方式填充整个二维、三维或更高维的空间。

[1] 中的插图,展示第1、2、3阶的希尔伯特曲线

本文描述用 C 语言绘画这种曲线的方法和技巧。

(题图是一个三维的希尔伯特曲线雕塑,来自 http://mathbun.com/


1. 递归、坐标变换

最初,我在 [2] 找到一种基于矢量的绘画方法。它的原理是,一阶的希尔伯特曲线是一个「⊓」的形状,我们可以生成四个顶点,把顶点连接起来。而第二阶则是用第一阶的四个顶点位置,加上旋转和缩小,去画4 个「⊓」形状。只要把最后生成的顶点连接起来,就可以生成整个曲线:

左图是一阶,起点在左下角。右图黑线部分是4个「⊓」形状,它们的中心为左图的顶点位置,加上旋转 -90、0、0、90度、以及缩小为一半。最后把所有顶点连起来,增加了灰线部分。

我们用一个影像缓冲去绘画这些垂直、水平线,最后用 极简的 PNG 编码函数 svpng() 来存储。

#include "svpng.inc"
#define S 512

unsigned char img[S * S * 3];
float px = 0.0f, py = 0.0f;

void setpixel(int x, int y) { /* ... */}
void lineto(float tx, float ty) {/* ... */}

void hilbert(float x, float y, float xi, float xj, float yi, float yj, int n) {
    if (n) {
        hilbert(x,                   y,                    yi / 2,  yj / 2,  xi / 2,  xj / 2, n - 1);
        hilbert(x + xi / 2,          y + xj / 2,           xi / 2,  xj / 2,  yi / 2,  yj / 2, n - 1);
        hilbert(x + xi / 2 + yi / 2, y + xj / 2 + yj / 2,  xi / 2,  xj / 2,  yi / 2,  yj / 2, n - 1);
        hilbert(x + xi / 2 + yi,     y + xj / 2 + yj,     -yi / 2, -yj / 2, -xi / 2, -xj / 2, n - 1);
    }
    else
        lineto(x + (xi + yi) / 2, y + (xj + yj) / 2);
}

int main() {
    hilbert(0.0f, 0.0f, 0.0f, S, S, 0.0f, 4);
    svpng(fopen("hilbert.png", "wb"), S, S,img, 0);
}

4阶的结果(为方便起见,起点位于左上角):

注:setpixel()lineto()较锁碎,不在此展示,完整代码位于https://github.com/miloyip/misc/blob/master/hilbert/hilbert.c


2. 方向编码

第一个方法要计算每个顶点的位置,计算量较大。另一种绘画希尔伯特曲线的方式,是考虑它的 L-system:

变数: L, R
常数: F, +, -
公理: L
规则:
L → − R F + L F L + F R −
R → − L F + R F R + F L −
F : 向前
- : 右转90°
+ : 左转90°

我们可以看到,两条规则其实是对称的,我们只需要把在迭归时,从 L 变成 R时把旋转方向反转。我参考了[3]的代码实现,但每次向前移动只绘画两个像素,当中p为笔的当前像素位置,参数d为当前方向,r记录旋转:

#include "svpng.inc"
#include <stdlib.h>

unsigned char* img, *p;
const int n = 4, s = (2 << n) - 1;

void step(int d) {
    int a[] = { 3, s * 3, -3, s * -3 }, i;
    for (i = 0; i < 2; i++, p += a[d & 3])
        p[0] = p[1] = p[2] = 255;
}

void hilbert(int d, int r, int n) {
    if (n--) {
        hilbert(d + r, -r, n); step(d + r);
        hilbert(d,      r, n); step(d    );
        hilbert(d,      r, n); step(d - r);
        hilbert(d - r, -r, n);
    }
}

int main() {
    p = img = calloc(s * s, 3);
    hilbert(0, 1, n);
    p[0] = p[1] = p[2] = 255;
    svpng(fopen("hilbert2.png", "wb"), s, s, img, 0);
    free(img);
}
放大 16 倍后的4阶结果

这种实现方式较简单,也能用最少的影像尺寸存储结果。


3. 简单文本

既然可以用较小的画布尺寸,我们直接把每个像素变换成** 两个字符:

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

char* img, *p;
const int n = 4, s = (2 << n) - 1, a[] = { 2, s * 2, -2, -s * 2 };

void step(int d) {
    int i;
    for (i = 0; i < 2; i++, p += a[d & 3])
        p[0] = p[1] = '*';
}

void hilbert(int d, int r, int n) { /* ... */ }

int main() {
    memset(p = img = malloc(s * s * 2), ' ', s * s * 2);
    hilbert(0, 1, n);
    p[0] = p[1] = '*';
    for (p = img; p < img + s * s * 2; p += s * 2)
        printf("%.*s\n", s * 2, p);
    free(img);
}

输出:


4. ASCII 字符美化

上面的输出有点丑。我在 [4] 里找到一种美化方法,可通过方向的改变来决定用什么字符,该映射为:

  • _
  • |_
  • |
  • |
  • __

例如:

我通过记录上一次的方向,与当前要步进的方向,去做这个映射。例如之前是向下的(l=3)当前要向右(d=0),即 的情况,便在笔的位置写进|_这两个字符。针对 4\times4=16 种组合,我用字符串表示:

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

char* img, *p;
const int n = 4 , s = 1 << n, a[] = { 2, s * 2, -2, -s * 2 };
int l = 3;

void step(int d) {
    d &= 3;
    p[0] = "_  ||||   _|   |"[l * 4 + d];
    p[1] = "_   _    ____   "[l * 4 + d];
    p += a[d];
    l = d;
}

void hilbert(int d, int r, int n) { /* ... */ }

int main() {
    p = img = malloc(s * s * 2);
    hilbert(0, 1, n);
    p[0] = p[1] = ' ';
    for (p = img; p < img + s * s * 2; p += s * 2)
        printf("%.*s\n", s * 2, p);
    free(img);
}

输出:

4阶
6阶

5. 无内存分配

为了尽一步简化,我希望能去掉画步的内存分配。方法是把逻辑改为,从上至下左至右打印字符,每画两个字符我们都遍历整个曲线,若曲线的位置和当前打印位置相同,才把那两个字符打印出来。这样会增加大量运算,但可以减少代码量。

#include <stdio.h>

const int n = 4, s = 1 << n, a[] = { 1, s, -1, -s };
int l = 3, p, q;

void step(int d) {
    d &= 3;
    if (p == q)
        printf("%.2s", &"__    | |_| |      ___|_ _    | "[l * 8 + d * 2]);
    p += a[l = d];
}

void hilbert(int d, int r, int n) { /* ... */ }

int main() {
    for (; q < s * s; q++, p = 0) {
        hilbert(0, 1, n);
        if (q % s == s - 1)
            putchar('\n');
    }
}

6. 代码压缩

最后,我把上面的代码压缩成 256 个字符:

const n=4,s=1<<n,a[]={1,s,-1,-s};
l=3,p,q;
t(d){d&=3;p-q||printf("%.2s",&"__    | |_| |      ___|_ _    | "[l*8+d*2]);p+=a[l=d];}
h(d,r,n){n--&&(h(d+r,-r,n),t(d+r),h(d,r,n),t(d),h(d,r,n),t(d-r),h(d-r,-r,n));}
main(){for(;p=0,q<s*s;++q%s||putchar(10))h(0,1,n);}
Try It Onlinetio.run

本文所有代码位于 https://github.com/miloyip/misc/tree/master/hilbert


相关文章

Milo Yip:用 C 语言画科赫雪花zhuanlan.zhihu.com图标Milo Yip:如何用 C 语言画这个图zhuanlan.zhihu.com图标如何用 C 语言画一个“圣诞树”?www.zhihu.com图标

参考

[1] D. Hilbert, Über die stetige Abbildung einer Linie auf ein Flächenstück, Mathematische Annalen 38 (1891), 459–460.

[2] http://www.soc.napier.ac.uk/~andrew/hilbert.html

[3] 形形色色的空间填充曲线 和 L-System

[4] ASCII Hilbert Curve