天泣記

2011-08-01 (Mon)

#1

しかし、ggplot2 だけでなく、plyr や reshape が気になる。

2011-08-06 (Sat)

#1

とりあえず、plyr と reshape の論文を読んでみた。

• Hadley Wickham, The Split-Apply-Combine Strategy for Data Analysis, Journal of Statistical Software, Vol. 40, Issue 1, Apr 2011
• Hadley Wickham, Reshaping Data with the reshape Package, Journal of Statistical Software, Vol. 21, Issue 12, Nov 2007

どちらもそれなりに面白かったが、どうも reshape の melt/cast が気になっていろいろと考えてみると、melt/cast はデータのレイヤ間の変換なのかもしれないと思えてきた。melt が低位レイヤへの変換で、cast が高位レイヤへの変換。

2011-08-08 (Mon)

#1

ふと思い立って、chkbuild が生成するログのファイルサイズをグラフにしてみる。

filesize.R:

library(ggplot2)
fs <- read.csv("2011-08/filesize.csv")
fs$time <- as.POSIXct(fs$base, "GMT", "%Y%m%dT%H%M%S")
print(qplot(time, size, data=fs, color=ext, size=I(1), ylim=c(0,1e6)))

ぬぅ。どんどん大きくなってやがる...

なお 2011-08/filesize.csv は chkbuild/tmp/public_html/ruby-trunk/log で ls -l した結果をてきとうにいじって作った以下のようなもの。

base,ext,size
20050718T131527,txt.gz,126211
20050718T134506,txt.gz,126234
20050718T144835,txt.gz,126423
...
20110808T060400Z,log.html.gz,519375
20110808T060400Z,log.txt.gz,517882
20110808T080500Z,diff.html.gz,915
20110808T080500Z,diff.txt.gz,549
20110808T080500Z,log.html.gz,519374
20110808T080500Z,log.txt.gz,517875
20110808T100500Z,diff.html.gz,1070
20110808T100500Z,diff.txt.gz,666
20110808T100500Z,log.html.gz,519721
20110808T100500Z,log.txt.gz,518166

#2

中身も計ってみるか。

chkbuild のログは "== name..." という行で各部分に分割できるので、分割してサイズを調べてみよう。

% ruby -rzlib -e '
def scan(fn)
Zlib::GzipReader.open(fn) {|f|
  base = fn.sub(/\..*/, "")           
  prevsec = nil
  prevpos = nil
  pos = 0
  f.each {|line|
    if /\A== ([a-z0-9]+)/ =~ line
      if prevsec
        puts [base, prevsec, pos - prevpos].join(",")
      end
      prevsec = $1
      prevpos = pos
    end
    pos += line.bytesize
  }
}
end
Dir["*.txt.gz"].grep(/[0-9Z]\.(log.txt.gz|txt.gz)/).sort.each {|arg|
  next if arg < "20060804"
  scan(arg)
}
'
20060804T123010,ruby,128
20060804T123010,start,37
20060804T123010,cvs,593008
20060804T123010,autoconf,51
20060804T123010,configure,10371
20060804T123010,make,41346
20060804T123010,version,88
20060804T123010,install,81887
20060804T123010,test,5571
20060804T123010,test,97762
20060804T123010,failure,56
20060804T131344,ruby,128
20060804T131344,start,37
20060804T131344,cvs,267227
20060804T131344,autoconf,51
20060804T131344,configure,10380
20060804T131344,make,49228
20060804T131344,version,88
20060804T131344,install,81772
20060804T131344,test,5582
20060804T131344,test,96349
20060804T131344,failure,56
20060804T140752,ruby,128
20060804T140752,start,37
20060804T140752,cvs,267550
20060804T140752,autoconf,51
...

base,section,size とヘッダをつけて、secsize.csv としておく。

secsize.R:

library(ggplot2)
secsize <- read.csv("2011-08/secsize.csv")
secsize$month <- substr(secsize$base, 1, 6)
secsize1 <- ddply(secsize, .(month), function(df)
        ddply(
          ddply(df, .(base, section), function(df2) data.frame(size=sum(df2$size))),
          .(section),
          function(df3) data.frame(size=mean(df3$size))))

print(qplot(month, size, data=secsize1, fill=section, color=I("black"),geom="bar") +
  scale_fill_manual(values=c("rubyspec"="red", "test"="blue", "svn"="yellow", "main"="cyan", "make"="magenta")) +
  scale_x_discrete(breaks=c("200701", "200801", "200901", "201001", "201101")))

テスト (rubyspec, test-all など) が大きいようだ。あと svn もけっこうある。

2011-08-15 (Mon)

#1

昨日までセキュリティ＆プログラミングキャンプに (特別チューターとして) 行ってきたのだが、その中でタイムトライアル処理系デバッグ演習 をやったときのメモを公開する。

ネタバレが嫌な人は見ないこと。

spcamp-timetrial.txt

問題文がないのでわかりにくいが、西尾さんの記事で触れられているものは (8) の部分である。test.rb というのは初耳で、与えられたものの中に含まれていたかどうか覚えていない。(演習の時は veryhard.zip だけが与えられて説明はとくになかった。)

ちなみに、各問の後に書いてある時刻は演習の残り時間で、veryhard.zip は 8分強で終わっている。いちおう演習時間中で (正確にはロスタイムに突入してたけど) 全部解いた。

#2

セキュリティ＆プログラミングキャンプの参加者に触発されて、Haskell の fixity について調べ、考えてみた。

(The Haskell 98 Report によれば) fixity は演算子の結合性と優先順位を宣言できる。たとえば infixr なら右結合である。

そして、この宣言は、変数の型宣言とかと同じような場所に書ける。つまり、変数のスコープと同じスコープで影響を与える。

とくに、where 節の中でも書けるのが剣呑である。

% cat a.hs
g a b c = a `f` b `f` c
  where
    f x y = x - y
    infixr `f`
% ghci a.hs
GHCi, version 6.8.2: http://www.haskell.org/ghc/  :? for help
Loading package base ... linking ... done.
[1 of 1] Compiling Main             ( a.hs, interpreted )
Ok, modules loaded: Main.
*Main> g 1 2 3
-4

`f` という二項演算子は、字面上、実際に使っている a `f` b `f` c よりも後に、右結合であることを宣言する infixr `f` が現れる。

つまり、最初から順に読んでいく限り、`f` を使っているところを読んだ時点でそれが右結合であることを知ることはできず、正しく構文木を構成するのは不可能である。

しかし、もし、宣言が事前に現れることが要求されるような仕様だったらどうだろう? where で後に宣言を書けるのは Haskell の (数学に由来する?) 悪癖だと思わないこともないので、そういう制限はありえると思う。そして、そういう仕様であれば、原理的には難しくないように思う。(あと、昔の A Gentle Introduction to Haskell には Fixity declarations must only appear at the very beginning of a Haskell module とあって、昔は実際そういう仕様だったようだ。)

yacc の類では、二項演算子の結合性や優先順位を宣言できるが、これは結局 shift/reduce conflict の解決のために使われる。逆にいえば、shift/reduce conflict の解決を自由に制御できれば、二項演算子の結合性や優先順位を変えられる。そのためには、shift/reduce conflict の解決を構文解析時まで遅延し、構文解析時に得られる情報を用いて解決の仕方を変えられるようにすればいい。

これは、原理としては、それほど難しくないのではないか。

実装としては、bison でそれを実現できるか、という点から始まり、できなければ拡張するか他の parser generator を探すか新しく作るかとか、まぁ、簡単にはいかないだろうけど。

2011-08-16 (Tue)

#1

サモアは今年末にタイムゾーンを 1日変えるらしい。

2011-12-29 の翌日が 2011-12-31 になるのだそうだ。

% TZ=Pacific/Samoa ruby -e 'p Time.now'         
2011-08-15 13:08:21 -1100

-1100 が +1300 になるということですかね。

まぁ、キリバスと同じような話なので目新しさはないな。

% TZ=Pacific/Kiritimati ruby -e 'p Time.local(1994,12,31,23,59,59)'
1994-12-31 23:59:59 -1000
% TZ=Pacific/Kiritimati ruby -e 'p Time.local(1994,12,31,23,59,59)+1'
1995-01-02 00:00:00 +1400

キリバスでは -1000 から +1400 だったか。1994-12-31 の翌日が 1995-01-02 になっている。

#2

W3C から Working with Time Zones なんて出てたのね。

2011-08-22 (Mon)

#1

iwlist コマンド (の scan サブコマンド) の結果をアクセスポイント毎に分割するのに Enumerable#slice_before を使った。

... String#split でも可能な気はする。

2011-08-28 (Sun)

#1

ふと正20面体を描こうと思い立って、OpenGL のサンプルをてきとうに改造したところこうなった。

rpol20.rb:

require "opengl"
require "glut"

PI = Math::PI
PI2 = Math::PI * 2

def ritri(a, b, c, ang1=:_, ang2=:_)
  # a**2 + b**2 = c**2
  # c*sin(ang1) = b
  # c*cos(ang1) = a
  # c*sin(ang2) = a
  # c*cos(ang2) = b
  ret = []
  if Numeric === a
    if Numeric === b
      ret << Math.hypot(a,b) if c == :*
      ret << Math.atan2(b,a)*360/PI2 if ang1 == :*
      ret << Math.atan2(a,b)*360/PI2 if ang2 == :*
    elsif Numeric === c
      ret << Math.sqrt(c*c - a*a) if b == :*
      ret << Math.acos(Float(a)/c)*360/PI2 if ang1 == :*
      ret << Math.asin(Float(a)/c)*360/PI2 if ang2 == :*
    elsif Numeric === ang1
      ret << a*Math.tan(ang1*PI2/360) if b == :*
      ret << a/Math.cos(ang1*PI2/360) if c == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << a/Math.tan(ang2*PI2/360) if b == :*
      ret << a/Math.sin(ang2*PI2/360) if c == :*
      ret << 90 - ang2 if ang1 == :*
    end
  elsif Numeric === b
    if Numeric === c
      ret << Math.sqrt(c*c - b*b) if a == :*
      ret << Math.asin(Float(b)/c)*360/PI2 if ang1 == :*
      ret << Math.acos(Float(b)/c)*360/PI2 if ang2 == :*
    elsif Numeric === ang1
      ret << b/Math.tan(ang1*PI2/360) if a == :*
      ret << b/Math.sin(ang1*PI2/360) if c == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << b*Math.tan(ang2*PI2/360) if a == :*
      ret << b/Math.cos(ang2*PI2/360) if c == :*
      ret << 90 - ang2 if ang1 == :*
    end
  elsif Numeric === c
    if Numeric === ang1
      ret << c*Math.cos(ang1*PI2/360) if a == :*
      ret << c*Math.sin(ang1*PI2/360) if b == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << c*Math.sin(ang2*PI2/360) if a == :*
      ret << c*Math.cos(ang2*PI2/360) if b == :*
      ret << 90 - ang2 if ang1 == :*
    end
  else
    given = [[:a, a], [:b, b], [:c, c], [:ang1, ang1], [:ang2, ang2]].reject {|var, val| !(Numeric === val) }
    given_vars = given.map {|var, val| var }
    raise ArgumentError, "information not enough (only #{given_vars.join(',')} given)"
  end
  ret
end

def myinit
  ambient = [0.0, 0.0, 0.0, 1.0]
  diffuse = [1.0, 1.0, 1.0, 1.0]
  position = [0.0, 3.0, 3.0, 0.0]

  lmodel_ambient = [0.2, 0.2, 0.2, 1.0]
  local_view = [0.0]

  GL::Light(GL::LIGHT0, GL::AMBIENT, ambient)
  GL::Light(GL::LIGHT0, GL::DIFFUSE, diffuse)
  GL::Light(GL::LIGHT0, GL::POSITION, position)
  GL::LightModel(GL::LIGHT_MODEL_AMBIENT, lmodel_ambient)
  GL::LightModel(GL::LIGHT_MODEL_LOCAL_VIEWER, local_view)

  GL::FrontFace(GL::CW)
  GL::Enable(GL::LIGHTING)
  GL::Enable(GL::LIGHT0)
  GL::Enable(GL::AUTO_NORMAL)
  GL::Enable(GL::NORMALIZE)
  GL::Enable(GL::DEPTH_TEST)
  GL::DepthFunc(GL::LESS)
end

def tri(p0, p1, p2)
  x0,y0,z0 = p0
  x1,y1,z1 = p1
  x2,y2,z2 = p2

  vx1, vy1, vz1 = x1 - x0, y1 - y0, z1 - z0
  vx2, vy2, vz2 = x2 - x0, y2 - y0, z2 - z0

  nx = vy1 * vz2 - vz1 * vy2 
  ny = vz1 * vx2 - vx1 * vz2 
  nz = vx1 * vy2 - vy1 * vx2 
  n = Math.sqrt(nx*nx + ny*ny + nz*nz)
  nx /= n
  ny /= n
  nz /= n

  GL::Normal(nx, ny, nz)
  GL::Begin(GL::TRIANGLES)
  GL::Vertex(x0, y0, z0)
  GL::Vertex(x1, y1, z1)
  GL::Vertex(x2, y2, z2)
  GL::End()
end

def quad(p0, p1, p2, p3)
  tri(p0, p1, p2)
  tri(p2, p3, p0)
end

$rotmatrix = nil
def render(x, y, ambr, ambg, ambb, difr, difg, difb, specr, specg, specb, shine)
  GL::PushMatrix()
  GL::Translate(x, y, 0.0)
  GL::MultMatrix($rotmatrix) if $rotmatrix
  GL::Scale(7, 7, 7)

  GL::Material(GL::FRONT, GL::AMBIENT, [ambr, ambg, ambb, 1.0])
  GL::Material(GL::FRONT, GL::DIFFUSE, [difr, difg, difb, 1.0])
  GL::Material(GL::FRONT, GL::SPECULAR, [specr, specg, specb, 1.0])
  GL::Material(GL::FRONT, GL::SHININESS, shine * 128.0)

  a, b = ritri(0.5, :*, :*, :_, 36)
  h1, = ritri(b, :*, 1)
  c, = ritri(0.5, :*, 1)
  h2, = ritri(b-a, :*, c)
  top_y = h2/2 + h1
  top = [0, h1 + h2/2, 0]
  bot = [0, -h1 - h2/2, 0]
  5.times {|i|
    p1 = [b*Math.sin(i*PI2/5), h2/2, b*Math.cos(i*PI2/5)]
    p2 = [b*Math.sin((i+1)*PI2/5), h2/2, b*Math.cos((i+1)*PI2/5)]
    p3 = [b*Math.sin((i+0.5)*PI2/5), -h2/2, b*Math.cos((i+0.5)*PI2/5)]
    p4 = [b*Math.sin((i+1.5)*PI2/5), -h2/2, b*Math.cos((i+1.5)*PI2/5)]
    tri(top,p1,p2)
    tri(bot,p4,p3)
    tri(p2,p1,p3)
    tri(p3,p4,p2)
  }

  GL::PopMatrix()
end

display = proc {
  GL::ClearColor(1,1,1,0)
  GL::Clear(GL::COLOR_BUFFER_BIT | GL::DEPTH_BUFFER_BIT)
  render(8.0, 8.0,
         0.135, 0.2225, 0.1575,
         0.54, 0.89, 0.63,
         0.316228, 0.316228, 0.316228,
         0.1)
  GLUT.SwapBuffers
}

reshape = proc {|w, h|
  GL::Viewport(0, 0, w, h)
  GL::MatrixMode(GL::PROJECTION)
  GL::LoadIdentity()
  if (w <= h)
    GL::Ortho(0.0, 16.0, 0.0, 16.0 *  h.to_f / w,
      -10.0, 10.0)
  else
    GL::Ortho(0.0, 16.0 * w.to_f / h, 0.0, 16.0,
      -10.0, 10.0)
  end
  GL::MatrixMode(GL::MODELVIEW)
}

$prev_pos = nil

mouse = proc {|button, state, x, y|
  p [button, state, x, y]
  if state == GLUT::UP
    $prev_pos = nil
  end
}

motion = proc {|x, y|
  $prev_pos = [x,y] if !$prev_pos
  dx = x - $prev_pos[0]
  dy = y - $prev_pos[1]
  $prev_pos = [x,y]
  p [dx,dy]
  GL::PushMatrix()
  GL::LoadIdentity()
  GL::Rotate(Math.hypot(dx,dy), dy, dx, 0)
  GL::MultMatrix($rotmatrix) if $rotmatrix
  $rotmatrix = GL::GetDoublev(GL::MODELVIEW_MATRIX)
  GL::PopMatrix()
  GLUT.PostRedisplay
}

keyboard = proc {|key, x, y|
  case (key)
  when ?q, 27 # q, ESC
    exit(0)
  end
}

GLUT::Init()
GLUT::InitDisplayMode(GLUT::DOUBLE | GLUT::RGB | GLUT::DEPTH)
GLUT::CreateWindow()
myinit()
GLUT::ReshapeFunc(reshape)
GLUT::DisplayFunc(display)
GLUT.MouseFunc(mouse)
GLUT.MotionFunc(motion)
GLUT.KeyboardFunc(keyboard)
GLUT::MainLoop()

ふむ、OpenGL というのはこういうものだったのか。

2011-08-29 (Mon)

#1

ついでに正12面体。

rpol12.rb:

require "opengl"
require "glut"
require "pp"

PI = Math::PI
PI2 = Math::PI * 2

def ritri(a, b, c, ang1=:_, ang2=:_)
  # a**2 + b**2 = c**2
  # c*sin(ang1) = b
  # c*cos(ang1) = a
  # c*sin(ang2) = a
  # c*cos(ang2) = b
  ret = []
  if Numeric === a
    if Numeric === b
      ret << Math.hypot(a,b) if c == :*
      ret << Math.atan2(b,a)*360/PI2 if ang1 == :*
      ret << Math.atan2(a,b)*360/PI2 if ang2 == :*
    elsif Numeric === c
      ret << Math.sqrt(c*c - a*a) if b == :*
      ret << Math.acos(Float(a)/c)*360/PI2 if ang1 == :*
      ret << Math.asin(Float(a)/c)*360/PI2 if ang2 == :*
    elsif Numeric === ang1
      ret << a*Math.tan(ang1*PI2/360) if b == :*
      ret << a/Math.cos(ang1*PI2/360) if c == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << a/Math.tan(ang2*PI2/360) if b == :*
      ret << a/Math.sin(ang2*PI2/360) if c == :*
      ret << 90 - ang2 if ang1 == :*
    end
  elsif Numeric === b
    if Numeric === c
      ret << Math.sqrt(c*c - b*b) if a == :*
      ret << Math.asin(Float(b)/c)*360/PI2 if ang1 == :*
      ret << Math.acos(Float(b)/c)*360/PI2 if ang2 == :*
    elsif Numeric === ang1
      ret << b/Math.tan(ang1*PI2/360) if a == :*
      ret << b/Math.sin(ang1*PI2/360) if c == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << b*Math.tan(ang2*PI2/360) if a == :*
      ret << b/Math.cos(ang2*PI2/360) if c == :*
      ret << 90 - ang2 if ang1 == :*
    end
  elsif Numeric === c
    if Numeric === ang1
      ret << c*Math.cos(ang1*PI2/360) if a == :*
      ret << c*Math.sin(ang1*PI2/360) if b == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << c*Math.sin(ang2*PI2/360) if a == :*
      ret << c*Math.cos(ang2*PI2/360) if b == :*
      ret << 90 - ang2 if ang1 == :*
    end
  else
    given = [[:a, a], [:b, b], [:c, c], [:ang1, ang1], [:ang2, ang2]].reject {|var, val| !(Numeric === val) }
    given_vars = given.map {|var, val| var }
    raise ArgumentError, "information not enough (only #{given_vars.join(',')} given)"
  end
  ret
end

def myinit
  ambient = [0.0, 0.0, 0.0, 1.0]
  diffuse = [1.0, 1.0, 1.0, 1.0]
  position = [0.0, 3.0, 3.0, 0.0]

  lmodel_ambient = [0.2, 0.2, 0.2, 1.0]
  local_view = [0.0]

  GL::Light(GL::LIGHT0, GL::AMBIENT, ambient)
  GL::Light(GL::LIGHT0, GL::DIFFUSE, diffuse)
  GL::Light(GL::LIGHT0, GL::POSITION, position)
  GL::LightModel(GL::LIGHT_MODEL_AMBIENT, lmodel_ambient)
  GL::LightModel(GL::LIGHT_MODEL_LOCAL_VIEWER, local_view)

  GL::FrontFace(GL::CW)
  GL::Enable(GL::LIGHTING)
  GL::Enable(GL::LIGHT0)
  GL::Enable(GL::AUTO_NORMAL)
  GL::Enable(GL::NORMALIZE)
  GL::Enable(GL::DEPTH_TEST)
  GL::DepthFunc(GL::LESS)
end

def tri(p0, p1, p2)
  x0,y0,z0 = p0
  x1,y1,z1 = p1
  x2,y2,z2 = p2

  vx1, vy1, vz1 = x1 - x0, y1 - y0, z1 - z0
  vx2, vy2, vz2 = x2 - x0, y2 - y0, z2 - z0

  nx = vy1 * vz2 - vz1 * vy2 
  ny = vz1 * vx2 - vx1 * vz2 
  nz = vx1 * vy2 - vy1 * vx2 
  n = Math.sqrt(nx*nx + ny*ny + nz*nz)
  nx /= n
  ny /= n
  nz /= n

  GL::Normal(nx, ny, nz)
  GL::Begin(GL::TRIANGLES)
  GL::Vertex(x0, y0, z0)
  GL::Vertex(x1, y1, z1)
  GL::Vertex(x2, y2, z2)
  GL::End()
end

def quad(p0, p1, p2, p3)
  tri(p0, p1, p2)
  tri(p2, p3, p0)
end

def pent(p0, p1, p2, p3, p4)
  tri(p0, p1, p2)
  tri(p0, p2, p3)
  tri(p0, p3, p4)
end

$rotmatrix = nil
def render(x, y, ambr, ambg, ambb, difr, difg, difb, specr, specg, specb, shine)
  GL::PushMatrix()
  GL::Translate(x, y, 0.0)
  GL::MultMatrix($rotmatrix) if $rotmatrix
  GL::Scale(5, 5, 5)

  GL::Material(GL::FRONT, GL::AMBIENT, [ambr, ambg, ambb, 1.0])
  GL::Material(GL::FRONT, GL::DIFFUSE, [difr, difg, difb, 1.0])
  GL::Material(GL::FRONT, GL::SPECULAR, [specr, specg, specb, 1.0])
  GL::Material(GL::FRONT, GL::SHININESS, shine * 128.0)

  a, b = ritri(0.5, :*, :*, :_, 36)
  c, d = ritri(:*, :*, 1, 54)

  r1 = b
  r2 = b*2*d

  h1, = ritri(r2-b, :*, 1)
  h2 = h1 * (a+b) / (a+b-c)
  hm = (h1+h2)/2

  bot = []
  top = []
  5.times {|i|
    j = i + 0.5
    k = i + 1.0
    l = i + 1.5
    t = [Math.sin(i*PI2/5), Math.cos(i*PI2/5)]
    u = [Math.sin(j*PI2/5), Math.cos(j*PI2/5)]
    v = [Math.sin(k*PI2/5), Math.cos(k*PI2/5)]
    w = [Math.sin(l*PI2/5), Math.cos(l*PI2/5)]

    pent([r1*w[0],hm,r1*w[1]],
         [r1*u[0],hm,r1*u[1]],
         [r2*u[0],hm-h1,r2*u[1]],
         [r2*v[0],hm-h2,r2*v[1]],
         [r2*w[0],hm-h1,r2*w[1]])

    pent([r1*t[0],-hm,r1*t[1]],
         [r1*v[0],-hm,r1*v[1]],
         [r2*v[0],-hm+h1,r2*v[1]],
         [r2*u[0],-hm+h2,r2*u[1]],
         [r2*t[0],-hm+h1,r2*t[1]])

    top << [r1*u[0], hm, r1*u[1]]
    bot << [r1*v[0], -hm, r1*v[1]]
  }

  pent(*top)
  pent(*bot.reverse)

  GL::PopMatrix()
end

display = proc {
  GL::ClearColor(1,1,1,0)
  GL::Clear(GL::COLOR_BUFFER_BIT | GL::DEPTH_BUFFER_BIT)
  render(8.0, 8.0,
         0.135, 0.2225, 0.1575,
         0.54, 0.89, 0.63,
         0.316228, 0.316228, 0.316228,
         0.1)
  GLUT.SwapBuffers
}

reshape = proc {|w, h|
  GL::Viewport(0, 0, w, h)
  GL::MatrixMode(GL::PROJECTION)
  GL::LoadIdentity()
  if (w <= h)
    GL::Ortho(0.0, 16.0, 0.0, 16.0 *  h.to_f / w,
      -10.0, 10.0)
  else
    GL::Ortho(0.0, 16.0 * w.to_f / h, 0.0, 16.0,
      -10.0, 10.0)
  end
  GL::MatrixMode(GL::MODELVIEW)
}

$prev_pos = nil

mouse = proc {|button, state, x, y|
  p [button, state, x, y]
  if state == GLUT::UP
    $prev_pos = nil
  end
}

motion = proc {|x, y|
  $prev_pos = [x,y] if !$prev_pos
  dx = x - $prev_pos[0]
  dy = y - $prev_pos[1]
  $prev_pos = [x,y]
  p [dx,dy]
  GL::PushMatrix()
  GL::LoadIdentity()
  GL::Rotate(Math.hypot(dx,dy), dy, dx, 0)
  GL::MultMatrix($rotmatrix) if $rotmatrix
  $rotmatrix = GL::GetDoublev(GL::MODELVIEW_MATRIX)
  GL::PopMatrix()
  GLUT.PostRedisplay
}

keyboard = proc {|key, x, y|
  case (key)
  when ?q, 27 # q, ESC
    exit(0)
  end
}

GLUT::Init()
GLUT::InitDisplayMode(GLUT::DOUBLE | GLUT::RGB | GLUT::DEPTH)
GLUT::CreateWindow()
myinit()
GLUT::ReshapeFunc(reshape)
GLUT::DisplayFunc(display)
GLUT.MouseFunc(mouse)
GLUT.MotionFunc(motion)
GLUT.KeyboardFunc(keyboard)
GLUT::MainLoop()

2011-08-30 (Tue)

#1

ひさしぶりに三次元を扱って、いろいろと検索していると、クォータニオン (四元数) が流行っているようである。(いや、流行っているというほどでもないか)

それなら利点があるのだろうと思ってしばらく検索しつつ考えると、結局、回転だけやりたければ 3x3 行列は自由度が高すぎるのがいかんのだろうという結論に達した。やりたくもないことができるのがよくないというのは理解できる。

しかし、回転だけであればクォータニオンでも自由度がひとつ多いのは許容範囲なのだろうか。

それはそれとして 3x3 行列からクォータニオンで失った自由度がなんなのかというのを考えると、方向によって異なる拡大縮小で 2, 剪断変形で (たぶん) 2 のような気がするのだが、クォータニオンの 4 とあわせて 8 にしかならない。3x3 行列は 9 なのでひとつ見逃しているはずだがなんだろう。もしかして剪断変形は 3 なのか?

... しばらく悩んでみたが、こう考えればいいのかもしれない。

3x3 行列を座標変換と考えると、変換元の3つの基底ベクトルに対応する変換先の3つのベクトル並べたものとみなせる。回転変換では、長さは変わらないので、変換先の3つのベクトルの長さはそれぞれ 1 である。また、角度も変わらないので、変換先の3つのベクトルは互いに直交する。というわけで、9個の値に、6個の制約があるので、回転の自由度は 3 である。

クォータニオンによる回転では、回転だけじゃなくて (方向性のない) 拡大縮小も可能なので、基底ベクトルの変換先の長さがそれぞれ 1 という制約が緩まり、同じ長さになるなら 1 でなくてもよいことになる。

というわけで、3x3 行列と比較してクォータニオンが失った自由度は、基底ベクトルの変換先の長さを異なる長さにする自由 (自由度 2) と、基底ベクトルを直交しないように変換する自由 (自由度 3) となる。

#2

クォータニオンを調べていると球面線形補間という話が出てくる。

spherical linear interpolation で slerp というようだ。

wikipedia (英語) の slerp に関する記述がわかりやすい。(クォータニオンとは関係なく定義できることが書いてある)

#3

ポインタのドラッグで多面体をぐるぐるまわせるようにしてあるのだが、それを実現するのは少し考える必要があった。(検索しても参考になるページをうまく見つけられなかった)

いまさら見つけたのだが、以下のページではそういうことを解説している。

<URL:http://marina.sys.wakayama-u.ac.jp/~tokoi/?date=20040321>

そのページではクォータニオンを使っているが、自分で書いたやつは回転行列を使っているという違いはある。でもドラッグの変位と直交する軸をとって回転するというやりかたは同じである。

2011-08-31 (Wed)

#1

で、正8面体。

rpol8.rb:

require "opengl"
require "glut"
require "pp"

PI = Math::PI
PI2 = Math::PI * 2

def ritri(a, b, c, ang1=:_, ang2=:_)
  # a**2 + b**2 = c**2
  # c*sin(ang1) = b
  # c*cos(ang1) = a
  # c*sin(ang2) = a
  # c*cos(ang2) = b
  ret = []
  if Numeric === a
    if Numeric === b
      ret << Math.hypot(a,b) if c == :*
      ret << Math.atan2(b,a)*360/PI2 if ang1 == :*
      ret << Math.atan2(a,b)*360/PI2 if ang2 == :*
    elsif Numeric === c
      ret << Math.sqrt(c*c - a*a) if b == :*
      ret << Math.acos(Float(a)/c)*360/PI2 if ang1 == :*
      ret << Math.asin(Float(a)/c)*360/PI2 if ang2 == :*
    elsif Numeric === ang1
      ret << a*Math.tan(ang1*PI2/360) if b == :*
      ret << a/Math.cos(ang1*PI2/360) if c == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << a/Math.tan(ang2*PI2/360) if b == :*
      ret << a/Math.sin(ang2*PI2/360) if c == :*
      ret << 90 - ang2 if ang1 == :*
    end
  elsif Numeric === b
    if Numeric === c
      ret << Math.sqrt(c*c - b*b) if a == :*
      ret << Math.asin(Float(b)/c)*360/PI2 if ang1 == :*
      ret << Math.acos(Float(b)/c)*360/PI2 if ang2 == :*
    elsif Numeric === ang1
      ret << b/Math.tan(ang1*PI2/360) if a == :*
      ret << b/Math.sin(ang1*PI2/360) if c == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << b*Math.tan(ang2*PI2/360) if a == :*
      ret << b/Math.cos(ang2*PI2/360) if c == :*
      ret << 90 - ang2 if ang1 == :*
    end
  elsif Numeric === c
    if Numeric === ang1
      ret << c*Math.cos(ang1*PI2/360) if a == :*
      ret << c*Math.sin(ang1*PI2/360) if b == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << c*Math.sin(ang2*PI2/360) if a == :*
      ret << c*Math.cos(ang2*PI2/360) if b == :*
      ret << 90 - ang2 if ang1 == :*
    end
  else
    given = [[:a, a], [:b, b], [:c, c], [:ang1, ang1], [:ang2, ang2]].reject {|var, val| !(Numeric === val) }
    given_vars = given.map {|var, val| var }
    raise ArgumentError, "information not enough (only #{given_vars.join(',')} given)"
  end
  ret
end

def myinit
  ambient = [0.0, 0.0, 0.0, 1.0]
  diffuse = [1.0, 1.0, 1.0, 1.0]
  position = [0.0, 3.0, 3.0, 0.0]

  lmodel_ambient = [0.2, 0.2, 0.2, 1.0]
  local_view = [0.0]

  GL::Light(GL::LIGHT0, GL::AMBIENT, ambient)
  GL::Light(GL::LIGHT0, GL::DIFFUSE, diffuse)
  GL::Light(GL::LIGHT0, GL::POSITION, position)
  GL::LightModel(GL::LIGHT_MODEL_AMBIENT, lmodel_ambient)
  GL::LightModel(GL::LIGHT_MODEL_LOCAL_VIEWER, local_view)

  GL::FrontFace(GL::CW)
  GL::Enable(GL::LIGHTING)
  GL::Enable(GL::LIGHT0)
  GL::Enable(GL::AUTO_NORMAL)
  GL::Enable(GL::NORMALIZE)
  GL::Enable(GL::DEPTH_TEST)
  GL::DepthFunc(GL::LESS)
end

def tri(p0, p1, p2)
  x0,y0,z0 = p0
  x1,y1,z1 = p1
  x2,y2,z2 = p2

  vx1, vy1, vz1 = x1 - x0, y1 - y0, z1 - z0
  vx2, vy2, vz2 = x2 - x0, y2 - y0, z2 - z0

  nx = vy1 * vz2 - vz1 * vy2 
  ny = vz1 * vx2 - vx1 * vz2 
  nz = vx1 * vy2 - vy1 * vx2 
  n = Math.sqrt(nx*nx + ny*ny + nz*nz)
  nx /= n
  ny /= n
  nz /= n

  GL::Normal(nx, ny, nz)
  GL::Begin(GL::TRIANGLES)
  GL::Vertex(x0, y0, z0)
  GL::Vertex(x1, y1, z1)
  GL::Vertex(x2, y2, z2)
  GL::End()
end

def quad(p0, p1, p2, p3)
  tri(p0, p1, p2)
  tri(p2, p3, p0)
end

def pent(p0, p1, p2, p3, p4)
  tri(p0, p1, p2)
  tri(p0, p2, p3)
  tri(p0, p3, p4)
end

$rotmatrix = nil
def render(x, y, ambr, ambg, ambb, difr, difg, difb, specr, specg, specb, shine)
  GL::PushMatrix()
  GL::Translate(x, y, 0.0)
  GL::MultMatrix($rotmatrix) if $rotmatrix
  GL::Scale(9, 9, 9)

  GL::Material(GL::FRONT, GL::AMBIENT, [ambr, ambg, ambb, 1.0])
  GL::Material(GL::FRONT, GL::DIFFUSE, [difr, difg, difb, 1.0])
  GL::Material(GL::FRONT, GL::SPECULAR, [specr, specg, specb, 1.0])
  GL::Material(GL::FRONT, GL::SHININESS, shine * 128.0)

  r, = ritri(0.5, 0.5, :*)
  h1, = ritri(r, :*, 1)

  top = [0,h1,0]
  bot = [0,-h1,0]
  4.times {|i|
    j = i+1
    p1 = [r*Math.sin(i*PI2/4), 0, r*Math.cos(i*PI2/4)]
    p2 = [r*Math.sin(j*PI2/4), 0, r*Math.cos(j*PI2/4)]
    tri(p1,p2,top)
    tri(p2,p1,bot)
  }

  GL::PopMatrix()
end

display = proc {
  GL::ClearColor(1,1,1,0)
  GL::Clear(GL::COLOR_BUFFER_BIT | GL::DEPTH_BUFFER_BIT)
  render(8.0, 8.0,
         0.135, 0.2225, 0.1575,
         0.54, 0.89, 0.63,
         0.316228, 0.316228, 0.316228,
         0.1)
  GLUT.SwapBuffers
}

reshape = proc {|w, h|
  GL::Viewport(0, 0, w, h)
  GL::MatrixMode(GL::PROJECTION)
  GL::LoadIdentity()
  if (w <= h)
    GL::Ortho(0.0, 16.0, 0.0, 16.0 *  h.to_f / w,
      -10.0, 10.0)
  else
    GL::Ortho(0.0, 16.0 * w.to_f / h, 0.0, 16.0,
      -10.0, 10.0)
  end
  GL::MatrixMode(GL::MODELVIEW)
}

$prev_pos = nil

mouse = proc {|button, state, x, y|
  p [button, state, x, y]
  if state == GLUT::UP
    $prev_pos = nil
  end
}

motion = proc {|x, y|
  $prev_pos = [x,y] if !$prev_pos
  dx = x - $prev_pos[0]
  dy = y - $prev_pos[1]
  $prev_pos = [x,y]
  p [dx,dy]
  GL::PushMatrix()
  GL::LoadIdentity()
  GL::Rotate(Math.hypot(dx,dy), dy, dx, 0)
  GL::MultMatrix($rotmatrix) if $rotmatrix
  $rotmatrix = GL::GetDoublev(GL::MODELVIEW_MATRIX)
  GL::PopMatrix()
  GLUT.PostRedisplay
}

keyboard = proc {|key, x, y|
  case (key)
  when ?q, 27 # q, ESC
    exit(0)
  end
}

GLUT::Init()
GLUT::InitDisplayMode(GLUT::DOUBLE | GLUT::RGB | GLUT::DEPTH)
GLUT::CreateWindow()
myinit()
GLUT::ReshapeFunc(reshape)
GLUT::DisplayFunc(display)
GLUT.MouseFunc(mouse)
GLUT.MotionFunc(motion)
GLUT.KeyboardFunc(keyboard)
GLUT::MainLoop()

#2

正6面体。

rpol6.rb:

require "opengl"
require "glut"
require "pp"

PI = Math::PI
PI2 = Math::PI * 2

def ritri(a, b, c, ang1=:_, ang2=:_)
  # a**2 + b**2 = c**2
  # c*sin(ang1) = b
  # c*cos(ang1) = a
  # c*sin(ang2) = a
  # c*cos(ang2) = b
  ret = []
  if Numeric === a
    if Numeric === b
      ret << Math.hypot(a,b) if c == :*
      ret << Math.atan2(b,a)*360/PI2 if ang1 == :*
      ret << Math.atan2(a,b)*360/PI2 if ang2 == :*
    elsif Numeric === c
      ret << Math.sqrt(c*c - a*a) if b == :*
      ret << Math.acos(Float(a)/c)*360/PI2 if ang1 == :*
      ret << Math.asin(Float(a)/c)*360/PI2 if ang2 == :*
    elsif Numeric === ang1
      ret << a*Math.tan(ang1*PI2/360) if b == :*
      ret << a/Math.cos(ang1*PI2/360) if c == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << a/Math.tan(ang2*PI2/360) if b == :*
      ret << a/Math.sin(ang2*PI2/360) if c == :*
      ret << 90 - ang2 if ang1 == :*
    end
  elsif Numeric === b
    if Numeric === c
      ret << Math.sqrt(c*c - b*b) if a == :*
      ret << Math.asin(Float(b)/c)*360/PI2 if ang1 == :*
      ret << Math.acos(Float(b)/c)*360/PI2 if ang2 == :*
    elsif Numeric === ang1
      ret << b/Math.tan(ang1*PI2/360) if a == :*
      ret << b/Math.sin(ang1*PI2/360) if c == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << b*Math.tan(ang2*PI2/360) if a == :*
      ret << b/Math.cos(ang2*PI2/360) if c == :*
      ret << 90 - ang2 if ang1 == :*
    end
  elsif Numeric === c
    if Numeric === ang1
      ret << c*Math.cos(ang1*PI2/360) if a == :*
      ret << c*Math.sin(ang1*PI2/360) if b == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << c*Math.sin(ang2*PI2/360) if a == :*
      ret << c*Math.cos(ang2*PI2/360) if b == :*
      ret << 90 - ang2 if ang1 == :*
    end
  else
    given = [[:a, a], [:b, b], [:c, c], [:ang1, ang1], [:ang2, ang2]].reject {|var, val| !(Numeric === val) }
    given_vars = given.map {|var, val| var }
    raise ArgumentError, "information not enough (only #{given_vars.join(',')} given)"
  end
  ret
end

def myinit
  ambient = [0.0, 0.0, 0.0, 1.0]
  diffuse = [1.0, 1.0, 1.0, 1.0]
  position = [0.0, 3.0, 3.0, 0.0]

  lmodel_ambient = [0.2, 0.2, 0.2, 1.0]
  local_view = [0.0]

  GL::Light(GL::LIGHT0, GL::AMBIENT, ambient)
  GL::Light(GL::LIGHT0, GL::DIFFUSE, diffuse)
  GL::Light(GL::LIGHT0, GL::POSITION, position)
  GL::LightModel(GL::LIGHT_MODEL_AMBIENT, lmodel_ambient)
  GL::LightModel(GL::LIGHT_MODEL_LOCAL_VIEWER, local_view)

  GL::FrontFace(GL::CW)
  GL::Enable(GL::LIGHTING)
  GL::Enable(GL::LIGHT0)
  GL::Enable(GL::AUTO_NORMAL)
  GL::Enable(GL::NORMALIZE)
  GL::Enable(GL::DEPTH_TEST)
  GL::DepthFunc(GL::LESS)
end

def tri(p0, p1, p2)
  x0,y0,z0 = p0
  x1,y1,z1 = p1
  x2,y2,z2 = p2

  vx1, vy1, vz1 = x1 - x0, y1 - y0, z1 - z0
  vx2, vy2, vz2 = x2 - x0, y2 - y0, z2 - z0

  nx = vy1 * vz2 - vz1 * vy2 
  ny = vz1 * vx2 - vx1 * vz2 
  nz = vx1 * vy2 - vy1 * vx2 
  n = Math.sqrt(nx*nx + ny*ny + nz*nz)
  nx /= n
  ny /= n
  nz /= n

  GL::Normal(nx, ny, nz)
  GL::Begin(GL::TRIANGLES)
  GL::Vertex(x0, y0, z0)
  GL::Vertex(x1, y1, z1)
  GL::Vertex(x2, y2, z2)
  GL::End()
end

def quad(p0, p1, p2, p3)
  tri(p0, p1, p2)
  tri(p2, p3, p0)
end

def pent(p0, p1, p2, p3, p4)
  tri(p0, p1, p2)
  tri(p0, p2, p3)
  tri(p0, p3, p4)
end

$rotmatrix = nil
def render(x, y, ambr, ambg, ambb, difr, difg, difb, specr, specg, specb, shine)
  GL::PushMatrix()
  GL::Translate(x, y, 0.0)
  GL::MultMatrix($rotmatrix) if $rotmatrix
  GL::Scale(9, 9, 9)

  GL::Material(GL::FRONT, GL::AMBIENT, [ambr, ambg, ambb, 1.0])
  GL::Material(GL::FRONT, GL::DIFFUSE, [difr, difg, difb, 1.0])
  GL::Material(GL::FRONT, GL::SPECULAR, [specr, specg, specb, 1.0])
  GL::Material(GL::FRONT, GL::SHININESS, shine * 128.0)

  r, = ritri(0.5, 0.5, :*)

  top = []
  bot = []
  4.times {|i|
    j = i+1
    p1 = [r*Math.sin(i*PI2/4), -0.5, r*Math.cos(i*PI2/4)]
    p2 = [r*Math.sin(j*PI2/4), -0.5, r*Math.cos(j*PI2/4)]
    p3 = [r*Math.sin(j*PI2/4), 0.5, r*Math.cos(j*PI2/4)]
    p4 = [r*Math.sin(i*PI2/4), 0.5, r*Math.cos(i*PI2/4)]
    quad(p1,p2,p3,p4)
    bot << p1
    top << p3
  }
  quad(*top)
  quad(*bot.reverse)

  GL::PopMatrix()
end

display = proc {
  GL::ClearColor(1,1,1,0)
  GL::Clear(GL::COLOR_BUFFER_BIT | GL::DEPTH_BUFFER_BIT)
  render(8.0, 8.0,
         0.135, 0.2225, 0.1575,
         0.54, 0.89, 0.63,
         0.316228, 0.316228, 0.316228,
         0.1)
  GLUT.SwapBuffers
}

reshape = proc {|w, h|
  GL::Viewport(0, 0, w, h)
  GL::MatrixMode(GL::PROJECTION)
  GL::LoadIdentity()
  if (w <= h)
    GL::Ortho(0.0, 16.0, 0.0, 16.0 *  h.to_f / w,
      -10.0, 10.0)
  else
    GL::Ortho(0.0, 16.0 * w.to_f / h, 0.0, 16.0,
      -10.0, 10.0)
  end
  GL::MatrixMode(GL::MODELVIEW)
}

$prev_pos = nil

mouse = proc {|button, state, x, y|
  p [button, state, x, y]
  if state == GLUT::UP
    $prev_pos = nil
  end
}

motion = proc {|x, y|
  $prev_pos = [x,y] if !$prev_pos
  dx = x - $prev_pos[0]
  dy = y - $prev_pos[1]
  $prev_pos = [x,y]
  p [dx,dy]
  GL::PushMatrix()
  GL::LoadIdentity()
  GL::Rotate(Math.hypot(dx,dy), dy, dx, 0)
  GL::MultMatrix($rotmatrix) if $rotmatrix
  $rotmatrix = GL::GetDoublev(GL::MODELVIEW_MATRIX)
  GL::PopMatrix()
  GLUT.PostRedisplay
}

keyboard = proc {|key, x, y|
  case (key)
  when ?q, 27 # q, ESC
    exit(0)
  end
}

GLUT::Init()
GLUT::InitDisplayMode(GLUT::DOUBLE | GLUT::RGB | GLUT::DEPTH)
GLUT::CreateWindow()
myinit()
GLUT::ReshapeFunc(reshape)
GLUT::DisplayFunc(display)
GLUT.MouseFunc(mouse)
GLUT.MotionFunc(motion)
GLUT.KeyboardFunc(keyboard)
GLUT::MainLoop()

#3

正4面体。

rpol4.rb:

require "opengl"
require "glut"
require "pp"

PI = Math::PI
PI2 = Math::PI * 2

def ritri(a, b, c, ang1=:_, ang2=:_)
  # a**2 + b**2 = c**2
  # c*sin(ang1) = b
  # c*cos(ang1) = a
  # c*sin(ang2) = a
  # c*cos(ang2) = b
  ret = []
  if Numeric === a
    if Numeric === b
      ret << Math.hypot(a,b) if c == :*
      ret << Math.atan2(b,a)*360/PI2 if ang1 == :*
      ret << Math.atan2(a,b)*360/PI2 if ang2 == :*
    elsif Numeric === c
      ret << Math.sqrt(c*c - a*a) if b == :*
      ret << Math.acos(Float(a)/c)*360/PI2 if ang1 == :*
      ret << Math.asin(Float(a)/c)*360/PI2 if ang2 == :*
    elsif Numeric === ang1
      ret << a*Math.tan(ang1*PI2/360) if b == :*
      ret << a/Math.cos(ang1*PI2/360) if c == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << a/Math.tan(ang2*PI2/360) if b == :*
      ret << a/Math.sin(ang2*PI2/360) if c == :*
      ret << 90 - ang2 if ang1 == :*
    end
  elsif Numeric === b
    if Numeric === c
      ret << Math.sqrt(c*c - b*b) if a == :*
      ret << Math.asin(Float(b)/c)*360/PI2 if ang1 == :*
      ret << Math.acos(Float(b)/c)*360/PI2 if ang2 == :*
    elsif Numeric === ang1
      ret << b/Math.tan(ang1*PI2/360) if a == :*
      ret << b/Math.sin(ang1*PI2/360) if c == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << b*Math.tan(ang2*PI2/360) if a == :*
      ret << b/Math.cos(ang2*PI2/360) if c == :*
      ret << 90 - ang2 if ang1 == :*
    end
  elsif Numeric === c
    if Numeric === ang1
      ret << c*Math.cos(ang1*PI2/360) if a == :*
      ret << c*Math.sin(ang1*PI2/360) if b == :*
      ret << 90 - ang1 if ang2 == :*
    elsif Numeric === ang2
      ret << c*Math.sin(ang2*PI2/360) if a == :*
      ret << c*Math.cos(ang2*PI2/360) if b == :*
      ret << 90 - ang2 if ang1 == :*
    end
  else
    given = [[:a, a], [:b, b], [:c, c], [:ang1, ang1], [:ang2, ang2]].reject {|var, val| !(Numeric === val) }
    given_vars = given.map {|var, val| var }
    raise ArgumentError, "information not enough (only #{given_vars.join(',')} given)"
  end
  ret
end

def myinit
  ambient = [0.0, 0.0, 0.0, 1.0]
  diffuse = [1.0, 1.0, 1.0, 1.0]
  position = [0.0, 3.0, 3.0, 0.0]

  lmodel_ambient = [0.2, 0.2, 0.2, 1.0]
  local_view = [0.0]

  GL::Light(GL::LIGHT0, GL::AMBIENT, ambient)
  GL::Light(GL::LIGHT0, GL::DIFFUSE, diffuse)
  GL::Light(GL::LIGHT0, GL::POSITION, position)
  GL::LightModel(GL::LIGHT_MODEL_AMBIENT, lmodel_ambient)
  GL::LightModel(GL::LIGHT_MODEL_LOCAL_VIEWER, local_view)

  GL::FrontFace(GL::CW)
  GL::Enable(GL::LIGHTING)
  GL::Enable(GL::LIGHT0)
  GL::Enable(GL::AUTO_NORMAL)
  GL::Enable(GL::NORMALIZE)
  GL::Enable(GL::DEPTH_TEST)
  GL::DepthFunc(GL::LESS)
end

def tri(p0, p1, p2)
  x0,y0,z0 = p0
  x1,y1,z1 = p1
  x2,y2,z2 = p2

  vx1, vy1, vz1 = x1 - x0, y1 - y0, z1 - z0
  vx2, vy2, vz2 = x2 - x0, y2 - y0, z2 - z0

  nx = vy1 * vz2 - vz1 * vy2 
  ny = vz1 * vx2 - vx1 * vz2 
  nz = vx1 * vy2 - vy1 * vx2 
  n = Math.sqrt(nx*nx + ny*ny + nz*nz)
  nx /= n
  ny /= n
  nz /= n

  GL::Normal(nx, ny, nz)
  GL::Begin(GL::TRIANGLES)
  GL::Vertex(x0, y0, z0)
  GL::Vertex(x1, y1, z1)
  GL::Vertex(x2, y2, z2)
  GL::End()
end

def quad(p0, p1, p2, p3)
  tri(p0, p1, p2)
  tri(p2, p3, p0)
end

def pent(p0, p1, p2, p3, p4)
  tri(p0, p1, p2)
  tri(p0, p2, p3)
  tri(p0, p3, p4)
end

$rotmatrix = nil
def render(x, y, ambr, ambg, ambb, difr, difg, difb, specr, specg, specb, shine)
  GL::PushMatrix()
  GL::Translate(x, y, 0.0)
  GL::MultMatrix($rotmatrix) if $rotmatrix
  GL::Scale(12, 12, 12)

  GL::Material(GL::FRONT, GL::AMBIENT, [ambr, ambg, ambb, 1.0])
  GL::Material(GL::FRONT, GL::DIFFUSE, [difr, difg, difb, 1.0])
  GL::Material(GL::FRONT, GL::SPECULAR, [specr, specg, specb, 1.0])
  GL::Material(GL::FRONT, GL::SHININESS, shine * 128.0)

  r, = ritri(0.5, :_, :*, 30)
  h, = ritri(r, :*, 1)
  ytop = h*3/4
  ybot = -h/4

  top = [0,ytop,0]
  bot = []
  3.times {|i|
    j = i+1
    p1 = [r*Math.sin(i*PI2/3), ybot, r*Math.cos(i*PI2/3)]
    p2 = [r*Math.sin(j*PI2/3), ybot, r*Math.cos(j*PI2/3)]
    tri(p1,p2,top)
    bot << p1
  }
  tri(*bot.reverse)

  GL::PopMatrix()
end

display = proc {
  GL::ClearColor(1,1,1,0)
  GL::Clear(GL::COLOR_BUFFER_BIT | GL::DEPTH_BUFFER_BIT)
  render(8.0, 8.0,
         0.135, 0.2225, 0.1575,
         0.54, 0.89, 0.63,
         0.316228, 0.316228, 0.316228,
         0.1)
  GLUT.SwapBuffers
}

reshape = proc {|w, h|
  GL::Viewport(0, 0, w, h)
  GL::MatrixMode(GL::PROJECTION)
  GL::LoadIdentity()
  if (w <= h)
    GL::Ortho(0.0, 16.0, 0.0, 16.0 *  h.to_f / w,
      -10.0, 10.0)
  else
    GL::Ortho(0.0, 16.0 * w.to_f / h, 0.0, 16.0,
      -10.0, 10.0)
  end
  GL::MatrixMode(GL::MODELVIEW)
}

$prev_pos = nil

mouse = proc {|button, state, x, y|
  p [button, state, x, y]
  if state == GLUT::UP
    $prev_pos = nil
  end
}

motion = proc {|x, y|
  $prev_pos = [x,y] if !$prev_pos
  dx = x - $prev_pos[0]
  dy = y - $prev_pos[1]
  $prev_pos = [x,y]
  p [dx,dy]
  GL::PushMatrix()
  GL::LoadIdentity()
  GL::Rotate(Math.hypot(dx,dy), dy, dx, 0)
  GL::MultMatrix($rotmatrix) if $rotmatrix
  $rotmatrix = GL::GetDoublev(GL::MODELVIEW_MATRIX)
  GL::PopMatrix()
  GLUT.PostRedisplay
}

keyboard = proc {|key, x, y|
  case (key)
  when ?q, 27 # q, ESC
    exit(0)
  end
}

GLUT::Init()
GLUT::InitDisplayMode(GLUT::DOUBLE | GLUT::RGB | GLUT::DEPTH)
GLUT::CreateWindow()
myinit()
GLUT::ReshapeFunc(reshape)
GLUT::DisplayFunc(display)
GLUT.MouseFunc(mouse)
GLUT.MotionFunc(motion)
GLUT.KeyboardFunc(keyboard)
GLUT::MainLoop()

#4

正多面体を描いて、というか記述して思うのは、プログラム内で正多面体の対称性を活用できていないということである。いちおう、一本の軸で対称なぶんはループでまとめたが、正多面体にはもっと多くの対称性があり、それらを利用できていない。

対称性をもっと積極的に利用し、記述量を減らすような仕掛けはデザインできないものかな。

まぁ正多面体に限ると 5種類しかないので、3ビットで指定できてしまうわけだが、もうちょっと広い対象を扱えることを意図して。

... ワイソフ記号というのがそうなのか?

#5

いろいろなダイスを売っているところを見つけた。

3面体、 5面体、 7面体、 14面体、 16面体、 24面体、 30面体が気になる。

3面体は楕球の側面から平面を 3枚削り出した感じ。3角柱の両端を丸めたもの、とも表現できるか。

5面体は見にくいが、おそらく (面取りされた) 3角柱か?

7面体は5角柱。

14面体、16面体はよくある 10面体ダイスの面を増やした感じ。

24面体は... 立方体の各面に浅い四角錐を張り付けたような感じ。(なんで凧形24面体にしないのか、とは思う)

30面体は菱形30面体。

[latest]