前回、ヒルベルト曲線を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]
ヒルベルト曲線上の任意の点:
コードは以下:
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
コードは以下:
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()
結果画像:
水平走査線との比較:
結果画像:
浮動小数点の問題(余談):
ヒルベルト曲線による画像スキャン:
この画像の特徴:
- サイズ(縦×横×RGB):512×512×3=262144×3=786432
- ヒルベルト曲線オーダー数:9 (4^9=262144)
オーダー9によるフルスキャンのコード:
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()
ヒルベルト曲線を利用した画像縮小:
- オーダー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