自作CNCマシン・レーザーカッターについて: ヒルベルト曲線：Hilbert Curve（Python）

2022年2月18日金曜日

ヒルベルト曲線：Hilbert Curve（Python）

今回は、Pythonでヒルベルト曲線を描いてみます。

ヒルベルト曲線の性質上、再帰アルゴリズムで解く方法がありますが、今回は再帰アルゴリズムを使わず、そのまま通常ループで解いていこうと思います（個人的にはそのほうが分かりやすいので）。

追記：再帰アルゴリズムを後半に追加しました。

環境： Python 3.8.5、Jupyter notebook

上画像、ORDER=6の場合（4096個の点をつないだ線）。

ORDERを上げていけば、面を埋め尽くしてしまうという空間充填曲線の一つ。

また拡大縮小しても相似形となるフラクタルな性質も持ち合わせています。

アルゴリズム：

まずヒルベルト曲線は4点からなるコの字を基準として、それの組み合わせで描いていきます。

コの字の組み合わせとして、90回転が4パターン、そして時計回り、反時計回りがあるので、合計4x2=8パターンありますが、実際はそのうちの4パターンしか使いません。

今回は、座標(0,0)を開始点とし移動方向が上右下のコの字を基本形とし、そのほか右上左、左下右、下左上の4パターンを使います。

ORDER=1の場合、上画像のような4点になります。

ルール1：（X座標が右向き、Y座標が上向きの場合）

4点を(0, 1, 2, 3)の順番とする。
原点0：（0, 0）を左下にする。
原点(0, 0)から上(0, 1)、右(1, 1)、下(1, 0)という順番で描く。

＊尚、プログラミング環境によってY座標が下向きの場合は、原点が左上となり、下、右、上（上下反転）。

次にORDER=2の場合（赤破線はORDER=1）。

黒番号は全体の通し番号（0〜15）。

青番号はORDER=2におけるローカルな番号（ORDER=1の番号と同じ位置関係）。

1組目の青コの字（通し番号0〜3）は赤0のエリア

2組目の青コの字（通し番号4〜7）は赤1のエリア

3組目の青コの字（通し番号8〜11）は赤2のエリア

4組目の青コの字（通し番号12〜15）は赤3のエリア

ルール2：

点の数は4ⁿ個（コの字4点の組数：2ⁿ個）、n=ORDER
コの字（4点）ごとに（0, 1, 2, 3）の番号をつける。
前回のORDER=1（赤破線）を参照し、青（0,1,2,3）の順番を入れ替える。

赤0のエリア：　青(0,1,2,3)→青(0,3,2,1)、青1と青3を入れ替え
赤1のエリア：　青(0,1,2,3)→青(0,1,2,3)、入れ替えなし
赤2のエリア：　青(0,1,2,3)→青(0,1,2,3)、入れ替えなし
赤3のエリア：　青(0,1,2,3)→青(2,1,0,3)、青0と青2を入れ替え

前回のORDER=1（赤破線）を参照し、その順番で青コの字を配置する。

赤0:(0,0)のエリア内に青1組目を配置する。

赤(0, 0)に赤単位長2を掛けた値(0,0)を青の点に加算する（赤0の場合は結局0）。

赤1:(0,1)のエリア内に青2組目を配置する。

赤(0, 1)に赤単位長2を掛けた値(0,2)を青の点に加算する。

赤2:(1,1)のエリア内に青3組目を配置する。

赤(1, 1)に赤単位長2を掛けた値(2,2)を青の点に加算する。

赤3:(1,0)のエリア内に青4組目を配置する。

赤(1, 0)に赤単位長2を掛けた値(2,0)を青の点に加算する。

（例）通し番号7の場合：

赤1(0,1)のエリア内で青3(1,0)であるため

赤1(0,1)×赤単位長2 + 青(1,0)×青単位長1 = (0,1)×2 + (1,0)×1 = (1,2)

よって通し番号7は最終的に(1,2)の座標になる。

ORDER=Nの内容を得るには、ORDER=N-1を参照し、ORDER=N-1の内容を得るには、ORDER=N-2を参照するように遡っていくので、たしかに再帰アルゴリズムが向いていることが分かります。しかし、今回はこのまま続行。

上記手続きのコード：

N = 2**ORDER
X,Y = [], []                             
for i in range(N * N):                   
    U = [[0, 0], [0, 1], [1, 1], [1, 0]]
    px, py = 0.5, 0.5                   
    for j in reversed(range(ORDER)):    
        ID = i // 4**j % 4              
        L = 2**j                        
        if   ID == 0:
            U[1], U[3] = U[3], U[1]     
        elif ID == 3:
            U[0], U[2] = U[2], U[0]      
        px += U[ID][0] * L              
        py += U[ID][1] * L
    X.append(px)                        
    Y.append(py)

上記コードの説明（ORDER=3の場合）：

ORDER=3のため合計2³×2³=64点（0〜63）を順番に一つずつforループ処理させる。
基本となるコの字型（ORDER=1）の座標を、U=[[0,0],[0,1],[1,1],[1,0]]という配列にしておく。
ORDER=3のとき、縦横に8個ずつ点があるため単位長を1とすれば全体のサイズは8x8となり、2のORDER乗だけ大きくなる（最終的に正規化し1×1のサイズに変更可）。
px, pyは原点座標であり、最終的に最小単位長1の半分(0.5)のオフセットがあるため、最初にずらしておく。
ネストされたfor文でORDERの数だけループさせる。reversed()で反転させたループとなり、j=2,1,0となる。
ID=i//4**j%4は、ORDER=jによる各点の位置関係を0〜3の値で返す。ORDER=3の場合は4進数の3ビットとなる（以下に「4進数」の説明あり）。
L=2**jは、ORDER=jに応じた単位長。

j=2: L=4 (ORDER=1)
j=1: L=2 (ORDER=2)
j=0: L=1 (ORDER=3)

ID=0の時U[1]とU[3]、ID=3の時U[0]とU[2]の入れ替えてコの字の向きを変える。
px+=U[ID][0]*Lは、入れ替えによって得られたUの要素とORDER数に応じた単位長を掛け合わせた値が加算され、最終的な座標が得られる（pyについても同様）。

青がORDER=3の結果。緑がORDER=2、赤がORDER=1。

黒番号は、全体の通し番号（0〜63）。

どのORDERにおいても基本的に左下：0、左上：1、右上：2、右下：3（ORDER=1を基準とする）。

各番号に対する座標：0:(0,0)、1:(0,1)、2:(1,1)、3:(1,0)

例えば、通し番号45の場合（上から3、右から2）、

ORDER=1（赤）においては2のエリア(1,1)、ORDER=1の単位長L=4

ORDER=2（緑）においては3のエリア(1,0)、ORDER=2の単位長L=2

ORDER=3（青）においては1のエリア(0,1)、ORDER=3の単位長L=1

となるため、

（px,　py）= (1,1)×4 + (1,0)×2 + (0,1)×1 = (4,4) + (2,0) + (0,1) = (6,5)

となる。

同様に57の場合（右から4番目、下から2番目）の場合、

ORDER=1（赤）においては3のエリア(1,0)、ORDER=1の単位長L=4

ORDER=2（緑）においては0のエリア(0,0)、ORDER=2の単位長L=2

ORDER=3（青）においては1のエリア(0,1)、ORDER=3の単位長L=1

よって、

（px,　py）= (1,0)×4 + (0,0)×2 + (0,1)×1 = (4,0) + (0,0) + (0,1) = (4,1)

となる。

４進数：

仕組みとしては4進数であり、ORDERが１増えるごとに4進数の桁も一つ増えていくことになります。同時に桁が増えるごとに、その桁に対応した単位長も2ⁿになり、桁の値と単位長を掛け合わせて合算するとその地点の座標が求まります。

ORDER=3の場合は、IDは4進数による3ビット（以下の数値の並びを縦に見る）になり、

ORDER=1(j=2):ID=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1...]

ORDER=2(j=1):ID=[0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3,0,0,0,0,1,1,1,1...]

ORDER=3(j=0):ID=[0,1,2,3,0,1,2,3,0,1,2,3,0,1,2,3,0,1,2,3,0,1,2,3...]

ORDER=1は16個ごと、ORDER=2は4個ごと、ORDER=3は1個ごとに[0,1,2,3]の繰り返しになっています。ORDER=3の値の並びが最小単位のコの字に対応していますが、例えば13番目（0番目から数えて）の点は、縦に[0][3][1]（赤い列の部分）であり、

ORDER=1(j=2): ID=[0] 左下:(0,0)、単位長:L=2**j=4

ORDER=2(j=1): ID=[3] 右下:(1,0)、単位長:L=2**j=2

ORDER=3(j=0): ID=[1] 左上:(0,1)、単位長:L=2**j=1

となり、ORDER=1のスケールでは左下(0,0)、ORDER=2のスケールでは右下(1,0)、ORDER=3のスケールでは左上(0,1)ということであり、各スケールにおける区画から住所のように座標を割り出せます。

ただし、途中で基準となるコの字型U=[[0,0],[0,1],[1,1],[1,0]]の要素の入れ替えがあるので、あくまでIDの番号は上の並びであって、U[0]〜U[3]の中身の並びは少し変わります。

以下は入れ替え後のU[ID]の4進数3ビット。

ORDER=1(j=2):U[ID]=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1...]

ORDER=2(j=1):U[ID]=[0,0,0,0,3,3,3,3,2,2,2,2,1,1,1,1,0,0,0,0,1,1,1,1...]

ORDER=3(j=0):U[ID]=[0,1,2,3,0,3,2,1,0,3,2,1,2,3,0,1,0,3,2,1,0,1,2,3...]

ORDER=1は、[0,1,2,3]の順で各16個ずつ合計64個。

ORDER=2は、[0,3,2,1][0,1,2,3][0,1,2,3][2,1,0,3]の順で各4個ずつ合計64個。

ORDER=3は、[0,1,2,3][0,3,2,1][0,3,2,1][2,3,0,1][0,3,2,1][0,1,2,3]...という並び方になっています。

先ほどの13番目の点（赤い列の部分）は、U[1]とU[3]の入れ替えによって[0][1][3]に変わっています。

この座標を求めると、

ORDER=1(j=2):[0]:(0,0),単位長:L=2**j=4

ORDER=2(j=1):[1]:(0,1),単位長:L=2**j=2

ORDER=3(j=0):[3]:(1,0),単位長:L=2**j=1

ということになり、各座標に単位長を掛け合わせて合算すると、

(0,0)×4 + (0,1)×2 + (1,0)×1 = (0,0)+(0,2)+(1,0) = (1,2)

よって、13番目（0番目から数えて）の点の座標は、(x,y) = (1,2)になります。

この入れ替えを含めたパターンがヒルベルト曲線におけるアルゴリズムの特徴だと思います。フラクタルだけあってこれもまたあるパターンの繰り返しになっています。

再帰アルゴリズム（追加）：

せっかくなので再帰アルゴリムも試してみました（単にforループを外しただけです）。

Pythonのシステム上、再帰回数がデフォルト:1000回に設定されており、必要に応じて設定変更しなければいけないようです。

今回の場合：

オーダー4：1024回

オーダー5：4096回

オーダー6：16384回

オーダー7：65536回

オーダー8：262144回

オーダー9：1048576回

になります。

しかしながら、設定を挙げてもオーダー5にするとパンクしてしまいます。colab.googleでもオーダー5が限界だったので、Jupyter自体のメモリー設定を変えるか、メモ化などの工夫が必要かもしれません。面倒なので無理にこの方法を使わなくでもいいかもしれません。

ちなみにfunctools.lru_cacheのメモ化ツールを使う場合は、arrayやlistを含んだ関数はダメなので、それらを全部タプルに変換する必要があります（メモ化の場合、値そのものを辞書機能のキーに割り当てるため、その際にarrayやlistだとキーにならないためタプルに変換する必要がある）。そうやって試してみましたが、それでもダメでした。


import sys
sys.setrecursionlimit(10 ** 4)            # Increase recursive functions memory(default: 1000) 

def hilbert(i, L, U):
    if i == 4**Order:
        return P + 0.5                    # Final output when the iteration finishes, offset(x, y)=(0.5, 0.5)
    if L >= 0:
        ID = i // 4**L % 4                # Quarternary index number
        if   ID == 0:
            U[1], U[3] = U[3].copy(), U[1].copy() # Turn the U-shape in clockwise and steps in anti-clockwise
        elif ID == 3:
            U[0], U[2] = U[2].copy(), U[0].copy() # Turn the U-shape in anti-clockwise and steps in anti-clockwise
        P[:, i] += U[ID] * 2**L           # add each step to the array for the output
        return hilbert(i, L-1, U)         # go to the next L(L=Order-1, Order-2 ... 0, -1)
    
    U = array([[0,0],[0,1],[1,1],[1,0]])  # initialize U and L
    L = Ord - 1
    return hilbert(i+1, L, U)             # go to the next i

Order = 4                                 # Level of Order (fixed value)
i = 0                                     # The number of iteration(i=0, 1 ... 4**Order-1)
L = Ord - 1                               # Length of U-shape(L=Order-1, Order-2 ... 0, -1)
U = array([[0,0],[0,1],[1,1],[1,0]])      # Initial U-shape coodinates
P = zeros((2, 4**Order))                  # Initial empty array for the output

P = hilbert(i, L, U)                      # X, Y = P[0], P[1]

まず設定したOrderによって点の合計数（4**Order）が決まり、その回数だけ再帰ループすることを終了条件（最後に各点座標を含んだarray:Pにオフセット0.5を加えた値を返します）にします。

次に、iに対して参照する各オーダーの座標を積算していきます。その際も再帰ループでオーダーの分ループし、各オーダーにおける単位長Lをもとに座標を求めていきます。

このオーダーのループが終われば次の点（i+1）を代入した最後の再帰ループに移行します。

という二段階の再帰ループになっています。

i+1が終了条件に達すれば、そこで終了し、すべての点座標を含んだPを返します。

上記の二重再帰を少し変更してみました。以下のように最初の再帰をforループに変えてあります（再帰ループを使う意味が本末転倒）。


def hilbert(i):
    if i == 4**Order:
        return P + 0.5                    # Final output when the iteration finishes, offset(x, y)=(0.5, 0.5)
    U = array([[0,0],[0,1],[1,1],[1,0]])  # initialize U and L
    for j in reversed(range(Order)):      # loop for each Order-scale steps: j=Order-1, Order-2 ... 0 
        ID = i // 4**j % 4                # Quarternary index number
        if   ID == 0:
            U[1], U[3] = U[3].copy(), U[1].copy() # Turn the U-shape in clockwise and steps in anti-clockwise
        elif ID == 3:
            U[0], U[2] = U[2].copy(), U[0].copy() # Turn the U-shape in anti-clockwise and steps in anti-clockwise
        P[:, i] += U[ID] * 2**j           # add each step to the output array
    return hilbert(i+1)                   # go to the next i

Order = 5                                 # Level of Order (fixed value)
P = zeros((2, 4**Order))                  # Initial empty array for the output
i = 0                                     # The number of iteration(i=0, 1 ... 4**Order-1)
P = hilbert(i)                            # X, Y = P[0], P[1]

再帰が一つ減った分、オーダー5までは動くようになりましたが、それ以上だとパンクしてしまいます。徐々にもともとのforループだけのコードに近くなってきたので、無理に再帰を使わなくてもいいかもしれません。再帰なしのforループだけのコードならオーダー9でも動きます。

さらに以下は、全体をforループで回し、その内部でオーダーに応じたステップ数を計算する部分だけ再帰にしてみた場合です。この場合もオーダー5でパンクしてしまいます。


def hilbert(n, L, U):
    for i in range(n, 4**Order):
        if L >= 0:
            ID = i // 4**L % 4                        # Quarternary index number
            if   ID == 0:
                U[1], U[3] = U[3].copy(), U[1].copy() # Turn the U-shape in clockwise and steps in anti-clockwise
            elif ID == 3:
                U[0], U[2] = U[2].copy(), U[0].copy() # Turn the U-shape in anti-clockwise and steps in anti-clockwise
            P[:, i] += U[ID] * 2**L           # add each step to the array for the output
            return hilbert(n, L-1, U)         # go to the next L(next Order: 1, 2 ... Order-1)
        U = array([[0,0],[0,1],[1,1],[1,0]])  # initialize U and L
        L = Order - 1
        n += 1
    return P

Order = 4                                 # Level of Order (fixed value)
P = zeros((2, 4**Order))                  # Initial empty array for the output
n = 0                                     # The number of iteration(i=0, 1 ... 4**Order-1)
L = Order - 1                             # Length of U-shape(L=Order-1, Order-2 ... 0)
U = array([[0,0],[0,1],[1,1],[1,0]])      # Initial U-shape coodinates

P = hilbert(n, L, U)                      # X, Y = P[0], P[1]

結果、このアルゴリズムだと再帰は向いていないようです。とりあえず現状では、再帰なしの通常forループで行ったほうがよさげです。

全体的なコード：

アニメーション描画のコードも追加してあります。

自作CNCマシン・レーザーカッターについて