Pythonあれこれ（Vol.74掲載）

著者：飯尾 淳

本連載では「Pythonを昔から使っているものの、それほど使いこなしてはいない」という筆者が、いろいろな日常業務をPythonで処理することで、立派な「蛇使い」に育つことを目指します。その過程を温かく見守ってください。皆さんと共に勉強していきましょう。第4回は、フラクタル図形の一種である「ヒルベルト曲線」を描く方法を解説します。

シェルスクリプトマガジン Vol.74は以下のリンク先でご購入できます。

図1　1次のヒルベルト曲線を描くPythonコード

#!/usr/bin/env python

from turtle import *
SIZE = 300
p = [[-0.5, 0.5], [-0.5, -0.5], [0.5, -0.5], [0.5, 0.5]]
def screen_pos(x):
  return (x[0]*SIZE, x[1]*SIZE)
setup(width = 2*SIZE, height = 2*SIZE)
color('red')
width(5)
hideturtle()
for x in p:
  goto(screen_pos(x))
mainloop()

図4　1 次のヒルベルト曲線を描くPythonコードの改良版

#!/usr/bin/env python

from turtle import *
SIZE = 300
penup_flag = True
p = [[-0.5, 0.5], [-0.5, -0.5], [0.5, -0.5], [0.5, 0.5]]
def screen_pos(x):
  return (x[0]*SIZE, x[1]*SIZE)
setup(width = 2*SIZE, height = 2*SIZE)
color('red')
width(5)
hideturtle()
penup()
for x in p:
  goto(screen_pos(x))
  if penup_flag:
    pendown()
    penup_flag = False
onkey(exit, 'q')
listen()
mainloop()

図7　2 次のヒルベルト曲線を描くPythonコード

#!/usr/bin/env python

from turtle import *
SIZE = 300
penup_flag = True
tm = [
  [ [0.0, -0.5, -0.5], [-0.5, 0.0,  0.5], [0.0, 0.0, 1.0] ],
  [ [0.5,  0.0, -0.5], [ 0.0, 0.5, -0.5], [0.0, 0.0, 1.0] ],
  [ [0.5,  0.0,  0.5], [ 0.0, 0.5, -0.5], [0.0, 0.0, 1.0] ],
  [ [0.0,  0.5,  0.5], [ 0.5, 0.0,  0.5], [0.0, 0.0, 1.0] ]
]
e = [ [ 1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0] ]
p = [ [-0.5, 0.5], [-0.5, -0.5], [0.5, -0.5], [0.5, 0.5] ]
def affine_transform(m, p):
  return [ m[0][0] * p[0] + m[0][1] * p[1] + m[0][2],
           m[1][0] * p[0] + m[1][1] * p[1] + m[1][2] ]
def mat_mul(m0, m1):
  return [ [m0[0][0]*m1[0][0]+m0[0][1]*m1[1][0]+m0[0][2]*m1[2][0],
            m0[0][0]*m1[0][1]+m0[0][1]*m1[1][1]+m0[0][2]*m1[2][1],
            m0[0][0]*m1[0][2]+m0[0][1]*m1[1][2]+m0[0][2]*m1[2][2]],
           [m0[1][0]*m1[0][0]+m0[1][1]*m1[1][0]+m0[1][2]*m1[2][0],
            m0[1][0]*m1[0][1]+m0[1][1]*m1[1][1]+m0[1][2]*m1[2][1],
            m0[1][0]*m1[0][2]+m0[1][1]*m1[1][2]+m0[1][2]*m1[2][2]],
           [m0[2][0]*m1[0][0]+m0[2][1]*m1[1][0]+m0[2][2]*m1[2][0],
            m0[2][0]*m1[0][1]+m0[2][1]*m1[1][1]+m0[2][2]*m1[2][1],
            m0[2][0]*m1[0][2]+m0[2][1]*m1[1][2]+m0[2][2]*m1[2][2]] ]
def screen_pos(x):
  return (x[0]*SIZE, x[1]*SIZE)
setup(width = 2*SIZE, height = 2*SIZE)
color('red')
width(5)
hideturtle()
penup()
for m in tm:
  p2 = [ affine_transform(mat_mul(e, m), x) for x in p ]
  for x in p2: 
    goto(screen_pos(x))
    if penup_flag:
      pendown()
      penup_flag = False
onkey(exit, 'q')
listen()
mainloop()

図9　書き換えた2 次のヒルベルト曲線を描くPythonコード

#!/usr/bin/env python

from turtle import *
import numpy as np
SIZE = 300
penup_flag = True
tm = [ np.matrix(x) for x in [
  [ [0.0, -0.5, -0.5], [-0.5, 0.0,  0.5], [0.0, 0.0, 1.0] ],
  [ [0.5,  0.0, -0.5], [ 0.0, 0.5, -0.5], [0.0, 0.0, 1.0] ],
  [ [0.5,  0.0,  0.5], [ 0.0, 0.5, -0.5], [0.0, 0.0, 1.0] ],
  [ [0.0,  0.5,  0.5], [ 0.5, 0.0,  0.5], [0.0, 0.0, 1.0] ] ]
]
e = np.eye(3) # 3次の単位行列
p = [ np.matrix(x).T for x in [
      [-0.5,  0.5, 1.0], [-0.5, -0.5, 1.0],
      [ 0.5, -0.5, 1.0], [ 0.5,  0.5, 1.0] ] ]
def screen_pos(x):
  return (x[0,0]*SIZE, x[1,0]*SIZE)
setup(width = 2*SIZE, height = 2*SIZE)
color('red', 'yellow')
width(5)
hideturtle()
penup()
for m in tm:
  for x in [ e * m * x for x in p ]:
    goto(screen_pos(x))
    if penup_flag:
      pendown()
      penup_flag = False
onkey(exit, 'q')
listen()
mainloop()

図10　図9 の赤枠部分を置き換えるPythonコード

def hilbert(n, m):
  if n > 1:
    for m_tm in tm: hilbert(n-1, m * m_tm)
  else:
    for x in [ m * x.T for x in p ]:
      goto(screen_pos(x))
      global penup_flag
      if penup_flag:
        pendown()
        penup_flag = False
hilbert(3, e)

図13　1 ～ 7 次のヒルベルト曲線を重ねて描画するPythonコード

#!/usr/bin/env python

from turtle import *
import numpy as np
SIZE = 300
MARGIN = 50
penup_flag = True
settings = [
  { 'color': 'forestgreen', 'width': 5 },
  { 'color': 'navy',        'width': 4 },
  { 'color': 'purple',      'width': 3 },
  { 'color': 'brown',       'width': 2 },
  { 'color': 'red',         'width': 2 },
  { 'color': 'orange',      'width': 1 },
  { 'color': 'yellowgreen', 'width': 1 }
]
tm = [ np.matrix(x) for x in [
  [ [0.0, -0.5, -0.5], [-0.5, 0.0,  0.5], [0.0, 0.0, 1.0] ],
  [ [0.5,  0.0, -0.5], [ 0.0, 0.5, -0.5], [0.0, 0.0, 1.0] ],
  [ [0.5,  0.0,  0.5], [ 0.0, 0.5, -0.5], [0.0, 0.0, 1.0] ],
  [ [0.0,  0.5,  0.5], [ 0.5, 0.0,  0.5], [0.0, 0.0, 1.0] ] ] ]
e = np.eye(3)
p = [ np.matrix(x).T for x in [
      [-0.5,  0.5, 1.0], [-0.5, -0.5, 1.0],
      [ 0.5, -0.5, 1.0], [ 0.5,  0.5, 1.0] ] ]
def draw_line_to(x):
  goto(x[0,0]*SIZE, x[1,0]*SIZE)
  global penup_flag
  if penup_flag:
    pendown()
    penup_flag = False
def reset():
  penup()
  global penup_flag
  penup_flag = True
def hilbert(n, m):
  if n > 1:
    for m_ in tm: hilbert(n-1, m * m_)
  else:
    for q in [ m * p_ for p_ in p ]: draw_line_to(q)
setup(width = 2*SIZE+MARGIN, height = 2*SIZE+MARGIN)
hideturtle()
for i in range(len(settings)):
  reset()
  color(settings[i]['color'], 'white')
  width(settings[i]['width'])
  hilbert(i+1, e)
onkey(exit, 'q')
listen()
mainloop()