grbl1.1+Arduino CNCシールドV3.5+bCNCを使用中。
BluetoothモジュールおよびbCNCのPendant機能でスマホからもワイヤレス操作可能。
その他、電子工作・プログラミング、機械学習などもやっています。
MacとUbuntuを使用。

CNCマシン全般について:
国内レーザー加工機と中国製レーザー加工機の比較
中国製レーザーダイオードについて
CNCミリングマシンとCNCルーターマシンいろいろ
その他:
利用例や付加機能など:
CNCルーター関係:



*CNCマシンの制作記録は2016/04/10〜の投稿に書いてあります。


ラベル アルゴリズム の投稿を表示しています。 すべての投稿を表示
ラベル アルゴリズム の投稿を表示しています。 すべての投稿を表示

2022年2月22日火曜日

ヒルベルト曲線/Hilbert curve(その2):曲線上の任意点と画像縮小

 前回、ヒルベルト曲線をPythonで描きました。今回はヒルベルト曲線のフラクタルな性質についての実験を行おうと思います。

環境: Python 3.8.5、Jupyter notebook


実験内容:

  • ヒルベルト曲線上に任意の点を複数打ち、オーダーを変えたときにそれらの点の位置関係は保たれるかどうか?
  • ヒルベルト曲線のルートで画像をスキャンし、そのスキャンデータをもとにオーダーを下げて画像縮小を行う
  • ヒルベルト曲線のルートで画像を描く(アニメーション)


ヒルベルト曲線のアルゴリズム:

前回との違い:(アルゴリズムの内容については前回参照)
  • 関数化
  • listをnumpy.arrayに変更(直接四則計算可能)
  • arrayの要素を入れ替える部分でcopy()を使用
    • U[1], U[3] = U[3].copy(), U[1].copy()
  • 正規化(norm)とオフセット(offset)を引数に与える
  • 出力されるX, Y=Pがヒルベルト曲線の各点の座標(array)

def hilbert(ORDER, norm=True, offset=0.5):            # argument: order=1, 2, 3 ... n, offset=(x:0.5, y:0.5)
    N = 2**ORDER                                      # the number of rows or columns
    P = zeros((2, N * N))                             # prepare an array for (X, Y) coodinates output
    for i in range(N * N):
        U = array([[0, 0], [0, 1], [1, 1], [1, 0]])   # four points under ORDER=1
        V = array([offset, offset])                   # offset: starting point
        for j in reversed(range(ORDER)):              # ORDER iteration: j=ORDER-1, ORDER-2 ... 0
            ID = i // 4**j % 4                        # quaternary index
            L = 2**j                                  # unit length: L=2**(ORDER-1), 2**(ORDER-2) ... 1 
            if   ID == 0:
                U[1], U[3] = U[3].copy(), U[1].copy() # swap 1st with 3rd when ID is 0
            elif ID == 3:
                U[0], U[2] = U[2].copy(), U[0].copy() # swap 0th with 2nd when ID is 3
            V += U[ID] * L                            # add each step under the ORDER
        P[:, i] = V                                   # preserve each points
    if norm:
        P /= N                                        # X, Y coodinates are nomalized as from 0 to 1
    return P                                          # returns X, Y = P[0], P[1]


ヒルベルト曲線上の任意の点:

ヒルベルト曲線は折れ曲がっていながらも連続した一本の線であるので、全体を一直線に引き延ばして全長を1(正規化)にしておきます。そこに任意の地点を設定するわけですが、例えば0.1の地点は全体の道のりの10%の位置になるというわけです。
線上の任意の点(0~1の値)がオーダーを変えたときにどのような座標になるのかを調べます。

コードは以下:


def Arbitrary_point(ap, X, Y): # ap: arbitrary point between 0 and 1
    tl = len(X) - 1
    n = int(ap * tl)
    m = ap * tl % 1
    P = [X[-1], Y[-1]]
    if n + m < tl:
        dx = X[n+1] - X[n]
        dy = Y[n+1] - Y[n]
        if dx == 0:
            P = [X[n], Y[n] + dy * m]
        elif dy == 0:
            P = [X[n] + dx * m, Y[n]]
    return P

この任意の点は、値によってはヒルベルト曲線上の二点間の間にも打つことができるようにしてあります。オーダー1のヒルベルト曲線上の0.1~0.9(0.1刻み)の地点は、以下の画像のように表現されます。9個の地点を新たな線でつないでいます。




この方法で、複数のオーダーによるヒルベルト曲線を重ね合わせて、任意の点(複数)を打って、それぞれどの程度違うのかを観察してみます。

コードは以下:


ORDERS = [3,4,5,6]            # list of Orders for Hilvert curves
AP = arange(0.1, 1, 0.1)      # list of arbitrary ratios on the lines: [0.1, 0.2 ... 0.9]

plt.figure(figsize=(10,10))
plt.axis([0,1,0,1])
cmap = plt.get_cmap('tab10')  # color table

X, Y = [], []
for i in ORDERS:
    X.append([])
    Y.append([])
    X[-1], Y[-1] = hilbert(i)
    plt.plot(X[-1], Y[-1], alpha=0.2, color=cmap(i))
    AX, AY = array([Arbitrary_point(ap, X[-1], Y[-1]) for ap in AP]).T
    plt.plot(AX, AY, color=cmap(i), alpha=0.7, label=str(i))
    plt.scatter(AX, AY, color=cmap(i))
    if i == ORDERS[0]:
        for j in range(len(AP)):
            plt.text(AX[j], AY[j], s=' ap='+str(round(AP[j], 3)))

plt.legend()
plt.show()


結果画像:



ここでは、オーダー3、4、5、6のヒルベルト曲線を重ね合わせています。そして任意の点は、0.1から0.9まで0.1刻みで9点打ってあります。
結果から分かることは、ほぼ同じような形をしており、フラクタルな性質が伺えます。オーダーが低い3のときはさすがにズレが大きいですが、オーダーが上がるほど誤差が少なくなってきます。


水平走査線との比較:

モニターなどで使われている走査線方式で同じことをやってみるとどうなるか試してみます。走査線は水平かつ平行に左から右、上から下へ走ることとします。

結果画像:



走査線の場合は、横幅のピクセル数が変わってしまうためズレが大きくなってしまいます。0.5の地点だけはきりのいい値なので似たような位置に集中しています。
走査線と比較すれば、ヒルベルト曲線は異なるスケールに対してもズレが少ないという傾向が見てとれます。


浮動小数点の問題(余談):

今回、小数の計算を行っている際に計算式と結果が一致しないことがあり、調べてみると浮動小数点の誤差問題ということらしいです。これはコンピュータ全般における問題でいくつか改善方法などあるようです。

例えば:

0.2 % 0.1 = 0.0

になりますが、

0.3 % 0.1 = 0.09999999999999998 

になってしまいます。
ためしに、以下のように0.1刻みで出力させると、

for i in arange(0, 1, 0.1):
    print(i)

出力結果:

0.0
0.1
0.2
0.30000000000000004
0.4
0.5
0.6000000000000001
0.7000000000000001
0.8
0.9

になり、どうやら0.3、0.6、0.7にこの問題が発生しているようです。
ちなみに、0~0.9までの値を0.1で割りその余りを100倍する計算:

for i in arange(0, 1, 0.1):
    print(i % 0.1 * 100)

答えはすべて0になるはずですが、出力は以下:

0.0
0.0
0.0
2.7755575615628914e-15
0.0
9.999999999999998
5.551115123125783e-15
2.7755575615628914e-15
0.0
9.999999999999998

0.3、0.6、0.7の出力は小さい誤差がありround()で有効桁数まで丸めてしまえば大丈夫かもしれませんが、問題は出力が9.999999999999998になっている部分。

0.5 % 0.1 * 100 = 9.999999999999998
0.9 % 0.1 * 100 = 9.999999999999998

本来0なのに、ここまでずれるとround()を使っても計算結果が大きくずれてしまいます。ちなみに最後に100を掛ける代わりに事前に0.5と0.1に100を掛けておけば問題ありません。

(0.5 * 100) % (0.1 * 100) = 0.0

ということで、厳密に計算するためのDecimalモジュールもありますが、連続で掛算や割算あるいはモデューロを使うときは注意が必要そうです。


ヒルベルト曲線による画像スキャン:

ヒルベルト曲線のルートに沿ってピクセル単位で画像をスキャンし、一本の連続した色データ(RGB)をつくり、各オーダーに応じた密度の違いを利用して画像縮小を行います。


おなじみのレナさんの画像を使って実験してみます。

この画像の特徴:

  • サイズ(縦×横×RGB):512×512×3=262144×3=786432
  • ヒルベルト曲線オーダー数:9 (4^9=262144)
通常、画像を読み込んでmatplotlibで表示させると原点が左上(Y軸が下向き)になることがありますが、その場合は:

plt.imshow(image, origin='lower')

のように'lower'にすると原点が左下になります。
また、画像が上下反対になっている場合は、numpy.flipud(image)や読み込んだarrayをimage[::-1, :, :]にすることでY軸方向のみの反転が可能になります。


オーダー9によるフルスキャンのコード:

このスキャンによって各ピクセルの色情報(RGB値)と位置(x,y)が得られます。読み取る順序はヒルベルト曲線のルート順になります。

DEFAULT_ORDER = int(np.log2(IMG_W))                   # DEFAULT_ORDER=9: 2**9=512
N = 2**DEFAULT_ORDER                                  # the number of rows or columns: 512

HX, HY = hilbert(DEFAULT_ORDER, norm=False, offset=0) # ORDER=9, normalization: turned off, offset=0
HX = HX.astype(int)                                   # dtype: float to int 
HY = HY.astype(int)

HC = np.zeros((N * N, 3), dtype=int) # prepare an empty 3-channel serial array for color pixel data
for i in range(N * N):
    HC[i, :] = IMG[HY[i], HX[i], :]  # scan color pixel data along with the Hilbert curve route

plt.figure(figsize=(15, 2))          # show a part of the serial color data, len(HC)=262144 pixels
start, end = 1000, 1050              # start and end point to sample the serial data 
plt.imshow([HC[start:end, :]], origin='lower', extent=(0,end-start,0,1))
plt.show()

3チャンネルのRGB情報を含んだ一本のデータの一部を出力してみると以下のようになります。


HCというarrayに格納してあり、サイズは(262144, 3)あります。
上の画像はその一部としてインデックス番号1000番から1049番までのピクセル。
各ピクセルの位置はヒルベルト曲線のルートに沿って(HX,HY)の変数に格納されています。


ヒルベルト曲線を利用した画像縮小:

上記のスキャンからHX, HY, HCのデータを得たので、それらをもとにオーダー数を下げたヒルベルト曲線で画像縮小を行います。

  • オーダー1:2×2ピクセル
  • オーダー2:4×4ピクセル
  • オーダー3:8×8ピクセル
  • オーダー4:16×16ピクセル
  • オーダー5:32×32ピクセル
  • オーダー6:64×64ピクセル
  • オーダー7:128×128ピクセル
  • オーダー8:256×256ピクセル
  • オーダー9:512×512ピクセル(元画像と同等)

縮小方法:

  • 各ピクセルの位置はオーダーに応じたヒルベルト曲線の点座標を用いる
  • 色情報については、オーダー9のフルスキャンの一直線データを以下の間隔でスキップして当てはめていく。近傍ピクセルの色調補正なし。
    • オーダー1:65536個ずつスキップ
    • オーダー2:16384個ずつスキップ
    • オーダー3:4096個ずつスキップ
    • オーダー4:1024個ずつスキップ
    • オーダー5:256個ずつスキップ
    • オーダー6:64個ずつスキップ
    • オーダー7:16個ずつスキップ
    • オーダー8:4個ずつスキップ
    • オーダー9:1個ずつ、スキップなし(計:262144ピクセル)

以下がコード:


ORDER = 6                                      # reduction to ORDER=6: 64x64=4096 pixels
N = 2**ORDER

HX, HY = hilbert(ORDER, norm=False, offset=0)  # obtain x, y coodinates of the Hilbert curve
HX = HX.astype(int)                            # convert to dtype as integer
HY = HY.astype(int)

HIM = zeros((N, N, IMG_CH), dtype=int)         # prepare an empty array for ORDER=6
for i in range(N * N):
    HIM[HY[i], HX[i], :] = HC[i * 4**(DEFAULT_ORDER - ORDER), :] # read the every 64th pixels from the serial data 

上記コードはオーダー6の場合。64個飛ばしに色情報を読み込んでいます。全体的なサイズは64×64で画像幅では1/8、面積的には1/64まで縮小したことになります。


結果の画像(オーダー9からオーダー6に縮小した場合):



オーダーを変えれば多少位置情報は異なりますが、オーダー6であればヒルベルト曲線のフラクタルな性質から一応認識できるレベルにあります。


水平走査線で同様に読み込ませた場合:




走査線の場合は、先ほどの実験通り画像が崩れてしまいます。縦と横の両方をスキップさせれば表現可能ですが、ヒルベルト曲線の場合は一つの情報をスキップさせるだけで済みます。



上画像はヒルベルト曲線におけるオーダー3(上の段:64ピクセル)から2(下の段:16ピクセル)へ縮小したときのデータの対応表です。ヒルベルト曲線では単に4個ずつスキップしていけばいいので簡単です。



そしてこちらは、水平走査線の場合です(同じ64ピクセルから16ピクセルへの縮小)。画像の横方向で2個ずつスキップし、さらに縦方向のために16スキップ(4個セットが16スキップ)しなければいけません。計算式は特に難しいことはないのですが、ヒルベルト曲線ほど単純ではありません。


ヒルベルト曲線スキャンのオーダーごとの画像:

画像上部にオーダー、サイズが書いてあります。
尚オーダー9が等倍(512×512ピクセル)で、それ以下のオーダー8からの画像になります。















認識可能なのはオーダー6以上でしょうか。オーダー5も画像そのものを小さく表示すればある程度は認識可能かもしれません。






オーダー6、5、4を極端に小さく表示してみました。
元々、何かのアイコン程度の画素数しかないので、このくらい小さく表示すればあまり違いが目立たなくなります。それでもオーダー4は厳しいかもしれません。



ヒルベルト曲線画像スキャンのアニメーション:

画像のピクセルを読み込むコードができたのでそのアニメーションもつくってみました。

オーダー6であれば4分程度でJupyter notebook上には表示できますが、それでもmatplotlibのデフォルトのメモリー設定を超えてしまうので、メモリー設定を変更しなければなりません。オーダー6のGIFアニメもつくりましたが、4096枚ものGIF画像を使用しているためかGIFファイルだけでも23MBあります。このブログにアップロードする際にエラーが出てしまったので、代わりにmp4の動画をアップロードしました。
また、オーダー7以上になると時間がかかりすぎるので試していません。
コード内容は以下のGistを参照してください。










関連:

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:

  • 点の数は4n個(コの字4点の組数:2n個)、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のため合計23×23=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進数:

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

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ループで行ったほうがよさげです。


全体的なコード:

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






関連:

2019年12月27日金曜日

クォータニオン / Quaternion:立方体の回転操作(その3)

前回までは単位球に対する一点の回転操作でしたが、今回はもう少し具体的にということで立方体の回転操作を試してみました(基本的なクォータニオン回転操作についてはこちらへ)。
立方体の8個の頂点を個別に計算して回転操作しています。
また、scale変数で立方体の縦横高さの比率を調整できるようにし、pos変数で位置の基準点を変えられるようにしています。

手順としては:
・8頂点(ベクトル)のそれぞれのスカラー値を求める
・8頂点(ベクトル)をスカラー値で正規化
・正規化した8頂点(ベクトル)を回転操作
・正規化されている8頂点(ベクトル)をスカラー倍してもとのスケールに戻す

ipywidgetsでインタラクティブに角度を変えられるようになっています。また立方体の縦横高さの比率も変更可。
オンライン上で実行するにはBinderで可能です(方法についてはこちらを参考に)。

以下がコード(Gist):
・回転角度と立方体の縦横高さ比が変更可能なバージョン
・回転角度と回転軸(単位ベクトル)が変更可能なバージョン
の二つがあります。

関連:
クォータニオン(四元数) / Quaternion / 回転制御(その1)
クォータニオン / 回転制御(その2):二点間から単位ベクトルと回転角度を求める


四元数
四元数
posted with amazlet at 19.12.22
今野 紀雄
森北出版
売り上げランキング: 67,351

2019年12月23日月曜日

クォータニオン / 回転制御(その2):二点間から単位ベクトルと回転角度を求める

前回に引き続きクォータニオン回転制御についてです。
今回は任意の二点(ベクトルAとB)だけが与られていて、AからBへ回転移動する際に必要な単位ベクトルと回転角度を求めてみます。

しかしながら、任意の二点(AとB)といっても、いくつかのケースがあります。

(1)AとBが単位球面上に存在する(|A|=|B|=1のとき)
(2)AとBのスカラー値は同じだけど、単位球面上に存在しない(|A|=|B|≠1のとき)
(3)AとBのスカラー値が異なり、単位球面上に存在しない(|A|≠|B|≠1)
4)そもそもAもBも単位ベクトルには無関係な任意の2点の場合
*(4)に関しては、AからBへの直線移動も含まれるし、単位ベクトル自体も任意に決定できるので今回は省略。このことから、任意の二点といってもあくまで単位球を設定しなければ逆算できないということがわかります。


(1)と(2)の場合(AとBのスカラー値が等しい場合):
(1)が一番簡単で、(2)も回転操作後に(1)の値を元のスカラー倍すればいいので、回転操作自体は同じ計算方法になります。以下はXZ平面において回転軸Yの場合。
単純にAが回転するとBになったという状況。回転操作のためには一旦ベクトルAを正規化してAnにしなければならないので、AnからBnへの回転計算をしたのち、正規化されているBnを元のスカラー倍してBとなります。

単位ベクトルを求める:
この場合、単位ベクトル(UV)は正規化されたAnとBnがあればクロス積(外積)で求めることができます(単位ベクトル(UV)は、ベクトルAnとベクトルBnに対して垂直/法線ベクトルとなるため)。ただし、Aのスカラー値とBのスカラー値は等しいという前提です(|A|=|B|)。
Pythonであれば、

A = [ax, ay, ax] # ベクトルA
B = [bx, by, bz] # ベクトルB

As = np.linalg.norm(A) # Aのスカラー値
Bs = np.linalg.norm(B) # Bのスカラー値 ただしAs=Bsという前提

An = np.array(A) / As  # 正規化
Bn = np.array(B) / Bs  # 正規化

UV = np.cross(An, Bn)  # 求められた単位ベクトル

回転角度を求める:
またAB間の回転角度(theta)はドット積(内積)とarccos()で求められ、

cos_AnBn = np.dot(An, Bn)
theta = np.arccos(cos_AnBn)

この二つの値(UV、theta)をクォータニオンの回転公式に代入し、Anを回転すれば計算結果がBnと一致します。そして正規化するまえの大きさに戻すため元のスカラー倍(Aのスカラー)します。

コード:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
# quaternion multiplication function
def QbyQ(q1, q2): # q = [w, x, y, z]
    w = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3]
    i = q1[0]*q2[1] + q1[1]*q2[0] + q1[2]*q2[3] - q1[3]*q2[2]
    j = q1[0]*q2[2] - q1[1]*q2[3] + q1[2]*q2[0] + q1[3]*q2[1]
    k = q1[0]*q2[3] + q1[1]*q2[2] - q1[2]*q2[1] + q1[3]*q2[0]
    return np.array([w, i, j, k])

# quaternion rotation function
def QPQc(P, UV, t): # P: quaternion, UV: unit vector, t: rotation angle theta
    q =  [np.cos(t/2),  UV[0] * np.sin(t/2),  UV[1] * np.sin(t/2),  UV[2] * np.sin(t/2)]
    qc = [np.cos(t/2), -UV[0] * np.sin(t/2), -UV[1] * np.sin(t/2), -UV[2] * np.sin(t/2)]
    qP = QbyQ(q, P)
    result = QbyQ(qP, qc)
    return result
# find unit vector and rotation angle for two vectors A and B
A = np.array([1, -1, 0.2])           # start point A
As = np.linalg.norm(A)               # scalar of A: |A|
An = A / As                          # normalization

B = np.array([0.4, -0.1, 1.1])       # end point B
Bs = np.linalg.norm(B)               # scalar of B: |B|
Bn = B / Bs                          # normalization

# find unit vector for An-Bn rotation
AnxBn = np.cross(An, Bn)             # orthogonal vector to An-Bn
UV = AnxBn / np.linalg.norm(AnxBn)   # unit vector (normalized): rotation axis 

# find rotation angle about An-Bn
# cos_t = (An dot Bb) / (|An|*|Bn|)  # (|An|*|Bn|) = An * Bn = 1 * 1 = 1
cos_AnBn = np.dot(An, Bn)
theta_AnBn = np.arccos(cos_AnBn)     # angle between An and Bn

# reproduct A-B rotation
Qa = np.insert(An, 0, 0)             # quaternion of An: [0, x, y, z]
Qb = QPQc(Qa, UV, theta_AnBn)        # quaternion rotation for An with UV and theta_AB
Bn_new = Qb[1:]                      # new Bn vector (normalized): [bx, by, bz]
B_new = Bn_new * Bs                  # new B vector

print('unit vector:', UV)
print('theta (degrees):', np.rad2deg(theta_AnBn))
print('initial vector B:', B)
print('new vector B    :', B_new) 


(3)の場合(AとBのスカラー値が異なる場合):
しかし(3)の場合は、単位球に対する任意の二点ABということから、Aを単純に回転させた結果がBになるとは限りません。以下のような状況。
Aを単純に回転させればA'(×印)にたどり着きますが、AとBのスカラー値(半径)が異なる場合はピンクの破線のような動きになります。Aは回転しつつ半径も縮まっていくという変化。
この場合、回転操作は(1)と同様に正規化されたベクトルAnで行い、最終的にAn(Bnと同じ)をBのスカラー倍(Bs倍)すればいいことになります。

回転軸となる単位ベクトルは(1)と同様にAnとBnのクロス積で求まります。

A = [ax, ay, ax] # ベクトルA
B = [bx, by, bz] # ベクトルB

As = np.linalg.norm(A) # Aのスカラー値
Bs = np.linalg.norm(B) # Bのスカラー値

An = np.array(A) / As
Bn = np.array(B) / Bs

UV = np.cross(An, Bn)  # 求められた単位ベクトル

AnとBnの回転角度に関しても(1)と同様にドット積とarccos()で求めます。

cos_AnBn = np.dot(An, Bn)
theta_AnBn = np.arccos(cos_AnBn)

求めた単位ベクトルUV(回転軸)とtheta_AnBn(回転角度)をクォータニオンの公式に代入すればBnが求まります。求めたBnに対して元のスカラー値Bsを掛ければBになります。
しかし、これでは移動前と移動後の二つの場面しかないので、途中の軌跡を描いて滑らかに移動しているかどうか確かめてみる必要がありそうです。つまり、先程の図のAからBへの移り変わり(ピンクの破線)を描いてみるというわけです。


AからBへの移動の中間補間:
まず移動のステップ数を決め、ステップ数に応じた回転角度とAからBへのスカラー値の変化量を決めます。

steps = 33
theta = np.linspace(0, theta_AnBn, steps)
ABs = np.lispace(As, Bs, steps)

これで、回転角度thetaは33ステップで0からtheta_AnBnまで変化します。同様にAB間のスカラー値ABsの変化も33ステップでAsからBsへ変化します。あとは公式に代入して33ステップ分の結果を得ます。

arw, arx, ary, arz = np.array([QPQc(Qa, UV, t) for t in theta]).T

上式で求まるのはAn-Bn間の回転結果(角度は違うけれどもスカラー値はどれも1)なので、各ステップに応じたスカラー値の変化量ABsを掛け合わせます。

ABx = arx * ABs
ABy = ary * ABs
ABz = arz * ABs

これで、各ステップごとに正規化されたベクトルに対して段階的にスカラー倍することになります。


ipywidgetsでインタラクティブ表示:
今回もipywidgetsでインタラクティブなグラフを表示してみます(ブログ上ではインタラクティブに動きませんが、mybinder.orgでこちらのgist URLを入力すればオンライン上でも動きます)。
関連:Binder:オンライン上でのJupyter Notebook(ipynbファイル)の実行
右端initial vector A(青)がinitialized vector B(緑)へ向かって移動する状態(ピンク)。unit vector:単位ベクトル(赤)が回転軸。
正規化されたAn-Bn間の回転移動は単位球面上にあります。AからBへの移動(ピンク破線)は単位球の同心円とずれているのがわかります。

クォータニオンの回転では、比較的簡単なパラメーターの調整だけで回転操作できることがわかりました。特にこのような移動途中の軌跡を描く際には便利なツールだと思います。

コード:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
# quaternion multiplication function
def QbyQ(q1, q2): # q = [w, x, y, z]
    w = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3]
    i = q1[0]*q2[1] + q1[1]*q2[0] + q1[2]*q2[3] - q1[3]*q2[2]
    j = q1[0]*q2[2] - q1[1]*q2[3] + q1[2]*q2[0] + q1[3]*q2[1]
    k = q1[0]*q2[3] + q1[1]*q2[2] - q1[2]*q2[1] + q1[3]*q2[0]
    return np.array([w, i, j, k])

# quaternion rotation function
def QPQc(P, UV, t): # P: quaternion, UV: unit vector, t: rotation angle theta
    q =  [np.cos(t/2),  UV[0] * np.sin(t/2),  UV[1] * np.sin(t/2),  UV[2] * np.sin(t/2)]
    qc = [np.cos(t/2), -UV[0] * np.sin(t/2), -UV[1] * np.sin(t/2), -UV[2] * np.sin(t/2)]
    qP = QbyQ(q, P)
    result = QbyQ(qP, qc)
    return result
def plot_base(elev=25, azim=-70):
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')
    ax.view_init(elev=elev, azim=azim)
    ax.set(xlim=(-1, 1), ylim=(-1 ,1), zlim=(-1, 1))
    ax.set(xlabel='X', ylabel='Y', zlabel='Z')
    ax.xaxis.pane.fill = ax.yaxis.pane.fill = ax.zaxis.pane.fill = False

    t = np.linspace(0, 2*np.pi, 128+1)
    alpha = 0.7
    ax.plot(np.cos(t), np.sin(t), [0]*len(t), linestyle=':', c='red', alpha=alpha)
    ax.plot(np.cos(t), [0]*len(t), np.sin(t), linestyle=':', c='red', alpha=alpha)
    ax.plot([0]*len(t), np.cos(t), np.sin(t), linestyle=':', c='red', alpha=alpha)
    ax.plot([-1, 1], [0, 0], [0, 0], linestyle=':', c='red', alpha=alpha)
    ax.plot([0, 0], [-1, 1], [0, 0], linestyle=':', c='red', alpha=alpha)
    ax.plot([0, 0], [0, 0], [-1, 1], linestyle=':', c='red', alpha=alpha)
    return ax
# find unit vector and rotation angle for two vectors A and B
A = np.array([1, -1, 0.2])           # start point A
As = np.linalg.norm(A)               # scalar of A: |A|
An = A / As                          # normalization

B = np.array([0.4, -0.1, 1.1])       # end point B
Bs = np.linalg.norm(B)               # scalar of B: |B|
Bn = B / Bs                          # normalization

# find unit vector for An-Bn rotation
AnxBn = np.cross(An, Bn)             # orthogonal vector to An-Bn
UV = AnxBn / np.linalg.norm(AnxBn)   # unit vector (normalized): rotation axis 

# find rotation angle about An-Bn
# cos_t = (An dot Bb) / (|An|*|Bn|)  # (|An|*|Bn|) = An * Bn = 1 * 1 = 1
cos_AnBn = np.dot(An, Bn)
theta_AnBn = np.arccos(cos_AnBn)     # angle between An and Bn

# reproduct A-B rotation
Qa = np.insert(An, 0, 0)             # quaternion of An: [0, x, y, z]
Qb = QPQc(Qa, UV, theta_AnBn)        # quaternion rotation for An with UV and theta_AB
Bn_new = Qb[1:]                      # new Bn vector (normalized): [bx, by, bz]
B_new = Bn_new * Bs                  # new B vector
# interactive plot with ipywidgets
from ipywidgets import interact

steps = 33
Arc = np.array([QPQc(Qa, UV, t) for t in np.linspace(0, theta_AnBn, steps)]).T # arc An-Bn
_, arx, ary, arz = Arc
ABs = np.linspace(As, Bs, steps) # scalar transition from A to B
ABx = arx * ABs                   # x term for arc A-B
ABy = ary * ABs                   # y term for arc A-B
ABz = arz * ABs                   # z term for arc A-B

@interact(idx=(0, steps-1, 1), elev=(0, 180, 1), azim=(-180, 180, 1))
def two_vec(idx=10, elev=25, azim=96):
    ax = plot_base(elev, -azim)

    ax.plot([0, A[0]], [0, A[1]], [0, A[2]], c='b', alpha=0.4)   # vector A
    ax.scatter(A[0], A[1], A[2], c='b', alpha=0.4)               # point A
    ax.text(A[0], A[1], A[2], s='  initial vector A')
    ax.plot([0, An[0]], [0, An[1]], [0, An[2]], c='b')           # vector An
    ax.scatter(An[0], An[1], An[2], c='b')                       # point An
    ax.text(An[0], An[1], An[2], s='normalized vector An  ', ha='right', va='top')

    ax.plot([0, UV[0]], [0, UV[1]], [0, UV[2]], c='r')           # unit vector
    ax.scatter(UV[0], UV[1], UV[2], c='r')
    ax.text(UV[0], UV[1], UV[2], s='unit vector  ', ha='right', va='top')

    ax.plot([0, B[0]], [0, B[1]], [0, B[2]], c='g', alpha=0.4)    # vector B
    ax.text(B[0], B[1], B[2], s='initialized vector B  ', ha='right', va='center')
    ax.scatter(B[0], B[1], B[2], c='g', alpha=0.4)                # vector B
    ax.plot([0, Bn[0]], [0, Bn[1]], [0, Bn[2]], c='g')            # vector Bn
    ax.scatter(Bn[0], Bn[1], Bn[2], c='g')                        # vector Bn
    ax.text(Bn[0], Bn[1], Bn[2], s='normalized vector Bn  ', ha='right', va='center')
    
    ax.plot([0, arx[idx]], [0, ary[idx]], [0, arz[idx]], color='magenta')           # vector An (move)
    ax.plot([0, arx[idx]*As], [0, ary[idx]*As], [0, arz[idx]*As], c='m', alpha=0.4) # move A    (move)
    ax.scatter(arx[idx]*As, ary[idx]*As, arz[idx]*As, c='m', alpha=0.4)             # point A   (move)
    ax.scatter(arx[idx], ary[idx], arz[idx], c='m')                                 # point An  (move)
    ax.plot(arx[:idx+1], ary[:idx+1], arz[:idx+1], linestyle='--', c='b', alpha=0.7)   # arc An-Bn
    ax.plot(arx[:idx+1]*As, ary[:idx+1]*As, arz[:idx+1]*As, ls='--', c='b', alpha=0.7) # arc A-A
    rate = As + (Bs - As) * idx / steps
    ax.scatter(arx[idx]*rate, ary[idx]*rate, arz[idx]*rate, c='m')                  # point A-B
    ax.plot(ABx[:idx+1], ABy[:idx+1], ABz[:idx+1], ls='--', c='m', alpha=0.9)       # arc A-B

関連:
クォータニオン(四元数) / Quaternion / 回転制御(その1)
クォータニオン / Quaternion:立方体の回転操作(その3)

四元数
四元数
posted with amazlet at 19.12.22
今野 紀雄
森北出版
売り上げランキング: 67,351

人気の投稿