Pygletでローレンツアトラクタ―を描く

Mark Summerfield氏の著書『実践Python3』でpygletを利用して3Dの描画を行うプログラムが紹介されていたので、Pygletでローレンツアトラクタ―を描くことを試みたので記録しておこう。
(本書の日本語版は斎藤康毅氏の翻訳で、私には「おう!あの『ゼロから作るDeep Learning』の著者か!」という第一印象が強かった。)
ローレンツアトラクタ―をC#で描画するプログラムを書いてみたが、透視投影は私には無理なので立体感がなく平板な印象のものになってしまった。
そこで本書で紹介されていたPyOpenGLかpygletを使用して3Dの描画を試みようと思った次第であるが、
PyOpenGLを使用したシリンダーを描画するスクリプトを実行すると、以下のようなエラーが出力された。

PS >python .\cylinder1.pyw
Traceback (most recent call last):
  File ".\cylinder1.pyw", line 139, in <module>
    main()
  File ".\cylinder1.pyw", line 22, in main
    glutInit(sys.argv)
  File "C:\Users\user1\AppData\Local\Programs\Python\Python37-32\lib\site-packages\OpenGL\GLUT\special.py", line 333, in glutInit
    _base_glutInit( ctypes.byref(count), holder )
  File "C:\Users\user1\AppData\Local\Programs\Python\Python37-32\lib\site-packages\OpenGL\platform\baseplatform.py", line 425, in __call__
    self.__name__, self.__name__,
OpenGL.error.NullFunctionError: Attempt to call an undefined function glutInit, check for bool(glutInit) before calling

glutInit関数が定義されていないということなのだが、確かにGLUTはインストールした記憶はない。
GLUTのインストールは選択肢は複数あり私にはちょっとハードルが高いなと感じた。一方、pygletの方はGLUTがなくても動くようだったのでpygletで試してみたら、立体感がある結果を得ることができた。

f:id:willwealth:20200122140611p:plain
Lorentz Attractor 3D by pyglet

計算部分がベタな感じのコードは以下の通り。

import os
import pyglet
from pyglet.gl import *
import numpy as np

PARAM_S = 10.0
PARAM_R = 50.0
PARAM_B = 8.0 / 3.0
DELTA_T = 0.001

START_X =  5.0
START_Y =  8.0
START_Z = 10.0
NUM_STEP  = 48000
INTERVAL = 10

SIZE = 800

ANGLE_INCREMENT = 1

COLOR = (255, 0, 0)
RADIUS = 0.5

def main():
    caption = "Lorenz Attractor( Sphere Version ) Drawed By Pyglet"
    width = height = SIZE
    resizable = True
    try:
        config = Config(sample_buffers=1, samples=4, depth_size=16,double_buffer=True)
        window = Window(width, height, caption=caption, config=config,resizable=resizable)
    except pyglet.window.NoSuchConfigException:
        window = Window(width, height, caption=caption,resizable=resizable)
    pyglet.app.run()

def vector(*args):
    return (GLfloat * len(args))(*args)

class Window(pyglet.window.Window):

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.set_minimum_size(SIZE, SIZE)
        self.xAngle = 2
        self.yAngle = -3
        self._initialize_gl()
        self._calc_attractor()

    def _initialize_gl(self):
        glClearColor(195/255, 248/255, 248/255, 1)
        glEnable(GL_DEPTH_TEST)
        glEnable(GL_POINT_SMOOTH)
        glHint(GL_POINT_SMOOTH_HINT, GL_NICEST)
        glEnable(GL_LINE_SMOOTH)
        glHint(GL_LINE_SMOOTH_HINT, GL_NICEST)
        glEnable(GL_COLOR_MATERIAL) # 1
        glEnable(GL_LIGHTING)
        glEnable(GL_LIGHT0)
        glLightfv(GL_LIGHT0, GL_POSITION, vector(0.5, 0.5, 1, 0))
        glLightfv(GL_LIGHT0, GL_SPECULAR, vector(0.5, 0.5, 1, 1))
        glLightfv(GL_LIGHT0, GL_DIFFUSE, vector(1, 1, 1, 1))
        glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 50)
        glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, vector(1, 1, 1, 1))
        glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE) # 2
        # 1 & 2 mean that we can use glColor*() to color materials

    def _draw_axes(self):
        glBegin(GL_LINES)                 # x-axis (traditional-style)
        glColor3f(1, 0, 0)
        glVertex3f( 0, 0, 0)
        glVertex3f( 1000, 0, 0)
        glEnd()

        glBegin(GL_LINES)                 # y-axis (traditional-style)
        glColor3f(0, 1, 0)
        glVertex3f(0, 0, 0)
        glVertex3f(0,  1000, 0)
        glEnd()

        glBegin(GL_LINES)                 # z-axis (traditional-style)
        glColor3f(0, 0, 1)
        glVertex3f(0, 0, 0)
        glVertex3f(0, 0,  1000)
        glEnd()

    def _calc_attractor(self):
        offset = np.array([ 0.19700702, 0.1629469,  45.41953137])
        self.pos_list = []
        x,y,z = START_X, START_Y, START_Z
        self.pos_list.append( np.array([x,y,z]) - offset )

        for i in range(1,NUM_STEP):
            xx = x + PARAM_S * ( -x + y ) * DELTA_T
            yy = y + ( -x * z + PARAM_R * x - y ) * DELTA_T
            zz = z + ( x * y - PARAM_B * z ) * DELTA_T
            if i % INTERVAL == 0 :
                self.pos_list.append( np.array([xx,yy,zz]) - offset )
            x,y,z = xx,yy,zz

    def on_draw(self):
        glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT)
        glMatrixMode(GL_MODELVIEW)
        glPushMatrix()
        glTranslatef(0, 0, -150)
        glRotatef(self.xAngle, 1, 0, 0)
        glRotatef(self.yAngle, 0, 1, 0)
        self._draw_axes()
        self._draw_spheres()
        glPopMatrix()
        pyglet.image.get_buffer_manager().get_color_buffer().save('screenshot.jpg')
        pyglet.image.get_buffer_manager().get_color_buffer().save('screenshot.png')

    def _draw_spheres(self):
        try:
            sphere = gluNewQuadric()
            gluQuadricNormals(sphere, GLU_SMOOTH)
            for x, y, z in self.pos_list:
                self._draw_sphere(sphere, x, y, z, RADIUS, COLOR)
        finally:
            gluDeleteQuadric(sphere)

    def _draw_sphere(self,sphere,x,y,z,radius,color):
        glPushMatrix()
        glTranslatef(x, y, z)
        glColor3ub(*color)
        gluSphere(sphere, radius, 24, 24)
        glPopMatrix()

    def on_resize(self, width, height):
        width = width if width else 1
        height = height if height else 1
        aspectRatio = width / height
        glViewport(0, 0, width, height)
        glMatrixMode(GL_PROJECTION)
        glLoadIdentity()
        gluPerspective(30.0, aspectRatio,  1.0, 1000.0)
        glMatrixMode(GL_MODELVIEW)
        glLoadIdentity()

# Escape is predefined to close the window so no on_key_press() needed
    def on_text_motion(self, motion): # Rotate about the x or y axis
        if motion == pyglet.window.key.MOTION_UP:
            self.xAngle -= ANGLE_INCREMENT
            #print(" xAngle = ",self.xAngle)
        elif motion == pyglet.window.key.MOTION_DOWN:
            self.xAngle += ANGLE_INCREMENT
            #print(" xAngle = ",self.xAngle)

        elif motion == pyglet.window.key.MOTION_LEFT:
            self.yAngle -= ANGLE_INCREMENT
            #print(" yAngle = ",self.yAngle)

        elif motion == pyglet.window.key.MOTION_RIGHT:
            self.yAngle += ANGLE_INCREMENT
            #print(" yAngle = ",self.yAngle)


if __name__ == "__main__":
    main()

# EOF

動画にするために少しづつ回転させた画像を作成するように変更を試みたが、同じスクリプトの中で実現することはできなかった。
そこでこのスクリプトを回転角度を引数で指定するように変更し、Powershellスクリプトで毎回、少しづつ回転させるという泥臭い方法で複数の画像を作成した。
しかし残念ながら、ffmpegで動画(mp4)を作成してみたが、ffplayでは再生できたがWindows Media Playerでは再生できなかった。
png形式の画像では駄目なのか、原因については不明のままであるが、もはや追究する気力もない。
次はUnityで試してみよう。

引数は辞書形式で与えて変数に反映。

    def analy_args(self):
        args = sys.argv[1:]
        param = dict([ i.split(':') for i in args ])
        print("param = {}".format(param))
        for key in ['rotx','roty']: # key of float value
            param[key] = float(param[key])
        self.xAngle = param['rotx']
        self.yAngle = param['roty']
        self.imgfile = param['img']

起動はPowerShellスクリプトで行い、Y軸回りに少しづつ回転させる(回転角度を大きくする)。

$prog = "python .\singleAttractor3D.py"
$count = 720
$dir = (Get-Date).ToString("yyyyMMdd")
$prefix = "img"

function execb ( $x, $y, $imgfile )
{
    "  x = " + $x
    "  y = " + $y
    "  imgfile = " + $imgfile
    $exec = $prog + " " + $x + " " + $y + " " + $imgfile
    $exec
    invoke-expression $exec
}

$time1 = get-date

for( $i = 0; $i -lt $count; $i ++)
{
    $f = "img:" + $prefix + $i.ToString("000#") + ".png"

    $a = "rotx:2";
    $v = 0.5 * $i
    $b = "roty:" + $v.ToString();

    execb $a $b $f
}

$time2 = get-date
$ts = new-timespan $time1 $time2
""
"Total Computation Time : {0:hh\:mm\:ss}" -f $ts

New-Item $dir -itemType Directory
Move-Item *.png $dir