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

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



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


ラベル Python の投稿を表示しています。 すべての投稿を表示
ラベル Python の投稿を表示しています。 すべての投稿を表示

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


全体的なコード:

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






関連:

2021年6月3日木曜日

IK(逆運動学):複数の円の交点から求める

 以前、CCD、FABRIK、ヤコビ行列などを用いて逆運動学を求めましたが、今回は複数の円を用いて、それらの交点から2D逆運動学を求めてみます。

環境:Python 3.8.5、Jupyter Notebook


これまでの方法:



今回の特長:

  1. 2D逆運動学
  2. 行列式、ヤコビ行列は使わない
  3. 繰り返し計算で目標値に近似する
  4. ガイドとなる円弧に沿ってアーム全体が配置される
  5. アームが交差しない
  6. ジョイント角が均等になる(各リンク長が等しい場合)
  7. 各リンク長が異なっても計算可能
  8. N個のリンクに対応

上図(アーム先端がターゲットに到達した状態):

  • 座標(0,0)がベース、LV[4]がアーム先端、TV(赤い×)がターゲット。
  • 青破線の円は各リンクの可動域、各ジョイント位置が円の中心。
  • 赤破線の円は各リンクを円弧状に配置するためのガイド、CPが円の中心。
  • 複数のリンクは赤破線の円に沿って配置されるため均等な角度となる(各リンク長が等しい場合)。


上図(初期状態:未到達時):計算手順

  • ベース(0,0)からターゲットTV(赤い×)まで線を引く(|TV|>0)
  • ベクトルTVの中点から垂直方向にCPを配置する(暫定的な配置)
  • (0,0)とTVを通り、CPを中心とした円を用意する(赤破線の円)、半径CR=|CP|
  • 各リンクのジョイント座標を赤破線の円上に配置する(青破線の円と赤破線の円の交点利用)
  • (0,0)、CP、TVの3点からなる角度Aを求める(下図)
  • (0,0)、CP、アーム先端(LV[4])の3点からなる角度Bを求める(下図)
  • 角度A+角度B=360度になるときアーム先端(LV[4])とTVが一致する
  • 角度A+角度Bが360度になるように、その差分Errorに応じてCPの座標を調節する計算式を用意する


Error = Thetas - tau
scaler += Error * 0.1
CP = np.array([-TV[1], TV[0]]) * scaler + TV/2

上記がCPの座標を調節する計算式。
Thetasは角度Aと角度Bの合計角、tauは2π=360度(目標値)、Errorはその差分。
CPの座標が変わることで赤破線の円の半径が変わる。常にCPはTV中点からの垂線上にある。
  • scaler=0のとき、TV/2がCPとなり、(0,0)とターゲット(TV)を直径とする円になる
  • scaler>0のとき、CPの位置はTV/2より上側になる
  • scaler<0のとき、CPの位置はTV/2より下側になる
つまり、目標値360度と角度A+角度Bの合計角の差分によってscaler値を増減させ、繰り返し計算によってアーム先端をターゲットに近似させていきます。

上図:

ターゲット(TV:赤い×)がベース(0,0)に近い場合でもアームが交差しない。ただし、|TV|>0の場合。


上図:

ターゲット(TV:赤い×)が遠く到達不可能な場合、ターゲット(TV)に向けてアームが伸びる。この場合、1000ループで計算終了。


上図:

リンク数10の場合。変数Nにリンク数を代入することで任意のリンク数に対応。


上図:

円の交点と角度の計算だけなので、各リンクが異なる長さの場合も計算可能。


コード:

  • Nはリンク数。TVはターゲットベクトル。LRは各リンク長(デフォルト:1.0)。
  • While文で繰り返し計算を行い、アーム先端(LV[4])とターゲットTVとの差分(Error)が0.0001以下、あるいは1000ループを超えるとループ解除。
  • intersects()関数は2つの円の中心座標と半径を引数にして、2つの交点を返しますが、今回は一方の交点(上側の交点座標)を利用。



2019年12月22日日曜日

クォータニオン(四元数) / Quaternion / 回転制御(その1)

今回はクォータニオンのアルゴリズムについてです(Python3.6/Jupyter notebook)。
クォータニオンは主に3DCG内のオブジェクトの回転制御で使われています。通常のXYZ軸ごとの角度制御よりも便利らしい。
しかしながら、二次元を複素平面で表現したように、三次元空間を四元数(一つの実数と三つの虚数i、j、k)によって表現しようとしている部分も興味深い。

クォータニオンは、

Q = w + x*i + y*j + z*k    (i2=j2=k2=ijk=-1)

で表現され、wは実数で残りの三つはXYZ軸に対応した複素数のようなもの。
Pythonにおいてはijkはないので、4つに分解して、

Q = [w, x, y, z]

という配列で表現しています。


回転制御:
通常回転制御する場合は、XYZ軸を何度ずつ動かせばいいのかと考えますが、クォータニオンの場合は回転軸一つと角度一つだけで制御します。

まず青いA点が単位球面上にあり、赤い回転軸(単位ベクトル:unit vector)を45度回すとBにたどり着くという状況があります。
この時点で、

・Aのベクトル(XYZ座標)
・回転軸(unit vector)のベクトル(XYZ座標)
・回転角(45度)

はすでに与えられており、この操作をするとBのベクトル(XYZ座標値)が得られるということになります。
言い換えれば、この操作方法によってAがどこに向かうのかを知りたいのであって、AからBにいくための操作方法を知りたいということではないので、きっかけが異なるとこの部分で少し戸惑います(AとBの二点間の回転移動を逆算して単位ベクトルと回転角度を求める内容は次回)。
しかしながら、赤い回転軸があれば単位球面上を最短距離でBに辿りつくことができます。わざわざ、X軸に何度、Z軸に何度、Y軸に何度と角度を三つに分解しなくても一つの角度だけで実現できるメリットがあり回転ツールとしては便利です。ただし3点はすべて単位曲面上にあるのでスカラー(半径)は1で座標値がそのままベクトルになります。

これをクォータニオンで計算するには、Aのベクトル[x, y, z]をクォータニオン式に変換し(wを付け足すだけ、w=0)、

A = [w, x, y, z]

とし、回転軸(unit vector)はベクトルのまま、

U = [ux, uy, uz]

として、クォータニオンの回転制御の公式である

w, x, y, z = Q*A*Q-

にいれれば、到着点であるBのベクトル/座標値[x, y, z]がわかります。
ただ、この公式の中身(QとQ-)が少し複雑です。

Q  = [cos(θ/2),  ux*sin(θ/2),  uy*sin(θ/2),  uy*sin(θ/2)]
Q- = [cos(θ/2), -ux*sin(θ/2), -uy*sin(θ/2), -uy*sin(θ/2)]

上式のθに回転角、回転軸:単位ベクトルのux, uy, uz値を代入。尚、Q-はQの共役クォータニオンで式中にマイナスがついており、Q*Q-=1という関係となっています。
さらに、クォータニオンの掛け算では各XYZにijkの虚数が含まれているため、実部wと虚部ijkの項に分けて計算します。

q1*q2 = [[q1w*q2w - q1x*q2x - q1y*q2y - q1z*q2z],
         [q1w*q2x + q1x*q2w + q1y*q2z - q1z*q2y],
         [q1w*q2y - q1x*q2z + q1y*q2w + q1z*q2x],
         [q1w*q2z + q1x*q2y - q1y*q2x + q1z*q2w]]
      = [w, x, y, z]

4段ありますが、上からw, x, y, z(実部, 虚部3つ)に対応しています。掛け算の結果として、[w, x, y, z]が得られるので、このうち[x, y, z]だけを取り出せばベクトル(XYZ座標値)になります。
また、ベクトル、内積、外積を使えばもう少し式が短くなります。

q1*q2:  w = q1w * q2w - (v1・v2)
        x, y, z = q1w * v2 + q2w * v1 + (v1×v2)

1段目が実部wに対応しており、2段目の虚部の式はXYZに対応した三つの値を返すので、そのままベクトルとして使えます。尚、(v1・v2)の部分はドット積(numpy.dot(v1, v2))、また(v1× v2)の部分はクロス積(numpy.cross(v1, v2))。ちなみに、v1=[q1x, q1y, q1z]はq1=[q1w, q1x, q1y, q1z]のXYZ成分です。
掛け算は複雑そうですが、コーディングする場合は一度関数にしておいて代入するだけなので楽です。
角度は一つしか登場してこないのでアニメーションなどつくるときには便利そうです。


ベクトルの正規化:
ベクトル[x, y, z]に実部wを加えればクォータニオンの式になりますが、先ほどの計算は全て正規化されていないとだめなので元のベクトルをスカラー値で割って正規化しておきます。

A = [x, y, z]

の場合、

√(x2 + y2 + z2) = 1

にならなければいけないので(ノルム:1)、

A_normalized = [x, y, z] / √(x2 + y2 + z2)

となります。実際Pythonではリストの割り算は一気にできないので、

A_normalized = np.array([x, y, z]) / ((x**2 + y**2 + z**2))**0.5
             = np.array(A) / np.linalg.norm(A)

になります。
同様にクォータニオンの場合は、

√(w2 + x2 + y2 + z2) = 1

でなければならないので、

Q = np.array([w, x, y, z]) / ((w**2 + x**2 + y**2 + z**2))**0.5

になります。
スカラー値で割って正規化して計算し、結果のベクトルを求め、それをスカラー倍すれば元のスケールに戻ります。


計算手順:
計算手順を再度まとめると、

A = [ax, ay, az] (これから動かそうとする点)

An = A / |A|     (正規化)

U = [ux, uy, uz] (単位ベクトル:回転軸)

Un = U / |U|     (正規化が必要であればしておく)

t = θ           (回転角)

以上が予め決めておくこと。
以下の式に角度と単位ベクトル(回転軸ベクトル)を代入しておきます。

Q  = [cos(t/2),  ux*sin(t/2),  uy*sin(t/2),  uy*sin(t/2)]
Q- = [cos(t/2), -ux*sin(t/2), -uy*sin(t/2), -uy*sin(t/2)]

そして、回転の公式で求めます。

w, x, y, z = Q*A*Q-  (QにAを掛けて、その結果にQ-を掛ける)

上式クォータニオン同志の掛け算は、以下の式を利用。

q1*q2 = [[q1w*q2w - q1x*q2x - q1y*q2y - q1z*q2z],
         [q1w*q2x + q1x*q2w + q1y*q2z - q1z*q2y],
         [q1w*q2y - q1x*q2z + q1y*q2w + q1z*q2x],
         [q1w*q2z + q1x*q2y - q1y*q2x + q1z*q2w]]
      = [w, x, y, z]

結果として[w, x, y, z]が求まるので、XYZ成分だけを取り出してベクトルにします。

[x, y, z]           (回転後のベクトル:単位球面上でのXYZ座標)

正規化される前のベクトルと同等のスカラー値に戻したい場合は、元のスカラー値を掛け直します。

[x, y, z] * |A|     (スカラー倍で元のスケールにする)

任意のベクトルA(青)を任意のunit vector(赤)を回転軸として45度動かしています。結果はベクトルB(緑)になります。

手順:
・まず、initial vector Aを正規化してnormalized vector Anにします。
・normalized vector Anをunit vector(回転軸)で45度回転させます。
・回転後はnormalized vector Bnが得られます。
・normalized vector BnにAのスカラー値をかけるとinitialized vector Bになります。
・よって、AからBに移動したことになります。

以下がPython(Jupyter Notebook)コード:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
A = np.array([1, -1, 1])             # initial vector A (arbitary vector)
As = np.linalg.norm(A)               # scalar of A: |A|
An = A / As                          # normalization: should be (sum(x**2 + y**2 + z**2))**0.5 = 1.0

U  = np.array([-0.1, -0.9, 0.2])     # unit vector (not normalized)
Us = np.linalg.norm(U)               # scalar of U: |U|
UV = U / Us                          # unit vector (normalized)

# 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

QA = np.insert(An, 0, 0)    # quaternion of An: [0, x, y, z]
theta = np.deg2rad(45)      # rotation angle
QB = QPQc(QA, UV, theta)    # quaternion rotation for An with UV and theta
Bn = QB[1:]                 # vectorize the quaternion result
B = Bn * As                 # new vector as same scalar as initial vector A
上記のプロット:
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

ax = plot_base()
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)
ax.text(A[0], A[1], A[2], s='  initial vector A', va='top')
ax.plot([0, An[0]], [0, An[1]], [0, An[2]], c='b')         # vector A normalized
ax.scatter(An[0], An[1], An[2], c='b')
ax.text(An[0], An[1], An[2], s='  normalized vector An', 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', 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='top')
ax.scatter(B[0], B[1], B[2], c='g', alpha=0.4)
ax.plot([0, Bn[0]], [0, Bn[1]], [0, Bn[2]], c='g')         # vector B normalized
ax.scatter(Bn[0], Bn[1], Bn[2], c='g')
ax.text(Bn[0], Bn[1], Bn[2], s='normalized vector Bn  ', ha='right', va='top')

Arc = [QPQc(QA, UV, t) for t in np.linspace(0, theta, 33)] # arc between A and B
arw, arx, ary, arz = np.array(Arc).T
ax.plot(arx, ary, arz, linestyle='--', c='b', alpha=0.6)
ax.plot(arx*As, ary*As, arz*As, linestyle='--', c='b', alpha=0.4)

まとめ:
クォータニオンの回転制御の公式に従えば簡単に操作できることがわかりました。それよりも3Dグラフに補助線を引いたりするコーディングのほうが大変でした。回転操作の計算内訳は少々複雑ですが、一旦ファンクション化してしまえば、あとは代入するだけなので手順はシンプルです。回転制御だけでなく、四元数そのものも気になります。
次回は、少し疑問に思っていた二点間の回転移動を逆算してみようと思います。

関連:
クォータニオン / 回転制御(その2):二点間から単位ベクトルと回転角度を求める
クォータニオン / Quaternion:立方体の回転操作(その3)

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

人気の投稿