前回、ヒルベルト曲線を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を参照してください。
関連：
