2018-06-07

Machinの公式を使ってπを近似する

数学小ネタ Ruby

どうも、drumathです。前回の記事からだいぶ日が経ってしまいましたね。
というのも、大学生活が始まり、なかなか執筆に時間がさけなかったのが原因です。
今日は気分で少し書いてみようと思ったのですが、テーマがまたπの近似という結果になってしまいました…（許して）

Machinの公式でπを求める！
終わりに

Machinの公式でπを求める！

冒頭で、大学生活が始まったと述べましたが、約1か月前の解析学の授業をきっかけにこのテーマを書こうと思いまして、
実はその授業の先生がπの近似値の世界記録を持ってるらしく（詳しい話はよく知りませんが...）、そんな経緯からなのでしょうか、逆三角関数の話をしているときにMachinの公式を教えてくださいました。
そもそも、逆三角関数は(辺の比)→(角度)という写像なので、πが絡む公式はいくつかあるらしく、下記のようなものが授業では紹介されました。

$\displaystyle \begin{align} 2 \tan^{-1} \frac{1}{2} - \tan^{-1} \frac{1}{7} = \frac{\pi}{4} \\ 2 \tan^{-1} \frac{1}{3} + \tan^{-1} \frac{1}{7} = \frac{\pi}{4} \\ 4 \tan^{-1} \frac{1}{5} - \tan^{-1} \frac{1}{239} = \frac{\pi}{4} \\ \end{align}$

実は、一番最後の「239！？」ってなるやつがMachinの公式です。

Machinの公式を示してみる

今からこの公式を使うわけですが、本当にこうなるのかとりあえず示します。
これには、

$\displaystyle \begin{align} \tan(\tan^{-1}x)=x \end{align}$

という性質を利用します。

まず、Machinの公式の両辺の正接をとります。すると

$\displaystyle \begin{align} \tan (4 \tan^{-1} \frac{1}{5} - \tan^{-1} \frac{1}{239} ) &= \tan \frac{\pi}{4} \\ &=1 \end{align}$

となりますね。よって、証明のゴールは左辺を変形していって１にすれば良いわけですね。やってみましょう！

$\displaystyle \begin{align} \tan (4 \tan^{-1} \frac{1}{5} - \tan^{-1} \frac{1}{239} ) について、\\ \tan^{-1} \frac{1}{5} = a, \tan^{-1} \frac{1}{239} = b とする。\\ \end{align}$

$\displaystyle \begin{align} また、\\ \tan 2a = \frac{2\tan a}{1-\tan^{2} a} &= \frac {2 \cdot \frac{1}{5} } { 1-(\frac {1} {5} )^{2} } \\ &= \frac{5}{12}\\ \\ \tan 4a = \frac{2\tan 2a}{1-\tan^{2} 2a} &= \frac{120}{119} \\ \\ よって\\ \tan(4a-b) = \frac{\tan 4a - \tan b}{1+ \tan 4a \tan b} &= \frac{ \frac{120}{119} - \frac{1}{239} } { 1 + \frac{120}{119} \cdot \frac{1}{239} } \\ &= 1 \end{align}$

とまぁ、最後にごりっと計算すると1になるわけです。すごいな。

πを計算するには

さて、証明ができたので、さっそく使っていきましょう
この公式を使ってどうやってπを求めるかというと、、、
積分だ！
待ってたぜぃって感じですが、どちらかというと積分は最後にちょこっとするだけで（というか積分ほぼやらないか…）、メインはマクローリン展開です。

まず、下の式を与えておきます。

$\displaystyle \begin{align} (\tan^{-1} x)^{'} = \frac{1}{1+x^{2}} \end{align}$

証明はめんどくさいのでやりません。ところで、彼(↑)はマクローリン展開すると

$\displaystyle \begin{align} \frac{1}{1+x^{2}} = 1-x^2+x^4-... \end{align}$

こうなりますよね。証明はしません。めんどくさいので。
これらがわかると下の式が導けますね。

$\displaystyle \begin{align} \tan^{-1} x &= \int (\tan^{-1} x)^{'} dx \\ &= \int \frac{1}{1+x^{2}} dx \\ &= \int (1-x^2+x^4-...)dx \\ &= x - \frac{x^3}{3} + \frac{x^5}{5} - ... \\ &= \lim_{n \to \infty} \sum^{n}_{k=0} (-1)^{k} \cdot \frac{x^{2k+1}}{2k+1} \end{align}$

おぉ、アークタンジェントがシグマで表せました！
ここで、

$\displaystyle \begin{align} f(x) = \lim_{n \to \infty} \sum^{n}_{k=0} (-1)^{k} \cdot \frac{x^{2k+1}}{2k+1} \end{align}$

とします。すると、Machinの公式より

$\displaystyle \begin{align} 4f \left( \frac{1}{5} \right) - f \left( \frac{1}{239} \right) = \frac{\pi}{4} \\ \\ \therefore \pi = 16f \left( \frac{1}{5} \right) - 4f \left( \frac{1}{239} \right) \end{align}$

プログラム書いてみた

さて、本記事も終盤です。実際にプログラミングで近似します。
その前に、プログラムのイメージがしやすいように(?)上でやったf(x)をちょっと改変したg(x,n)を以下のように定義します。

$\displaystyle \begin{align} g(x, n) = \sum^{n}_{k=0} (-1)^{k} \cdot \frac{x^{2k+1}}{2k+1} \end{align}$

あんまり見た目は変わりませんが、プログラムの構造がわかりやすいと思ったので…
gを使ってπを表すと、以下のようになります。

$\displaystyle \begin{align} \pi = \lim_{n \to \infty} 16g \left( \frac{1}{5}, n \right) - 4g \left( \frac{1}{239}, n \right) \end{align}$

となりますね。つまり、gという関数を実装して、二つ連結してループで回す、という感じです。

以下、私が書いたRubyのコードです。

def maclaurin(x, n)
  sum = 0.0
  0.upto(n) do |k|
    sum += (-1.0)**k * x**(2.0 * k + 1.0) / (2.0 * k + 1.0)
  end
  sum
end

def machin(n)
  4 * maclaurin(1.0 / 5.0, n) - maclaurin(1.0 / 239.0, n)
end

gets.chomp.to_i.times do |i|
  puts "#{i}: #{4 * machin(i)}"
end

ちなみに、rubocopからは「引数名が短い！」って怒られてますw

破壊力ぅ...

どうしてMachinの公式がπの近似に使われるか。
それはこの公式の破壊力がすさまじいからです(語彙力)
なんとサンプル数（ここでいうn）を10にすると、もうRubyではこれ以上の近似はできません！
収束がめっちゃ早い。

0: 3.18326359832636
1: 3.1405970293260603
2: 3.1416210293250346
3: 3.1415917721821773
4: 3.1415926824043994
5: 3.1415926526153086
6: 3.141592653623555
7: 3.1415926535886025
8: 3.141592653589836
9: 3.1415926535897922
10: 3.141592653589794
11: 3.141592653589794
12: 3.141592653589794
13: 3.141592653589794
14: 3.141592653589794
15: 3.141592653589794

終わりに

今回、英語の宿題を残して爆速で記事を書きましたが、久しぶりに楽しかったです。
またヒマになったら記事を上げていこうともうので、よかったらご覧になってください。
最後まで読んでいただき、ありがとうございました～。

2018-03-19

もはや自前鯖はいらない！Node.jsで作ったシンプルなチャットを開発～デプロイまでをブラウザだけでやる

今回は最近はまっているNode.jsについて書こうと思います。昔からリアルタイムチャットアプリを作るのは夢だったので、今回はそれをお題にしてみます。なお、同じ記事をqiitaにも投稿しました。

どんなものを、どうやって作るか

今回はこんなものを作ります。ルーム機能や、アカウント機能などない、シンプルなチャットです。ソケットの練習のために作ったみたいなものです。実際にデプロイしたものをこちらに公開しています。

どうやって作るの？

今回の開発では、リアルタイム通信機能が不可欠です。そのため、Socket.ioというパッケージを使います。Socket.ioと言ってしまったのでもうお分かりかもしれませんが、バックエンドはNode.jsを使います。オールJSで開発だ！いぇいということで、今回の開発プランはこうです。

Cloud9でBlankなワークスペースを作る
Node.jsで鯖を立てる
Socket.ioでソケット通信のイベントを記述（ここまでバックエンド）
クライアントサイドのhtmlを整える
クライアントサイドjsでソケットのイベント処理
Sassでスタイルを整える（ここまでフロントエンド）
Git周辺を整え、Herokuへデプロイ

ゆえに、今回はHTML、jsの基本はわかっている前提で、主にnode.jsやSocket.io、Herokuへのデプロイについて書きたいと思います。

Cloud9とは

Cloud9とは、ネット上で開発が行えるIDEで、Linux環境が使えます。私はWin機を使っているので、Linuxを使いたい時によく重宝しています。PythonもRubyもnodeもperlもphpも全部入ってるんでWebサービス開発にはもってこいの環境です。入門時は特に環境構築で手間取るので、入門者や、ちょっといじってみたいという方にお勧めです。過去にCloud9についての記事も書いたことがあるのでよろしければ読んでみてください。

開発だ！

Node.jsでHello

さて、まずはCloud9でワークスペースを「Blank」で作ります。すると、最初にREADMEが表示されると思いますが、別に消していただいていいでしょう。ではワークスペースにserver.jsを追加しましょう。フォルダのアイコンを右クリックか、ターミナルで

touch server.js

と打ってください。では、さっそくNode.jsで鯖を立てましょう。 server.jsに以下のようにコードを書きます。

let server = require("http").createServer((req, res)=>{
    res.writeHead(200, {'Content-Type': 'text/plane'})
    res.write('Hello')
    res.end()
})

server.listen(process.env.PORT, process.env.IP, ()=>{
    console.log("server")
})

ポート番号とホストはCloud9からprocess.env.PORTとprocess.env.IPを使えと言われます。そして上のほうにあるRunボタンを押すとnodeが走ってCloud9でWebページがホストされますので、見てみてください。Runボタンを押したときに教示されるターミナルの一番上に表示されているURLがそれです。「Hello」と表示されたでしょうか。

Expressの導入

それではフレームワークであるExpressとソケット通信のパッケージ、Socket.ioを導入しましょう。まず、Node.jsのパッケージを管理するnpmを初期化します。

npm init

質問されますので、その都度適当なものを入力します。するとpackage.jsonが作られます。そしたら

npm install --save express socket.io

とコマンドを打ってください。この２つがインストされるはずです。では、Expressを使ってHello,Worldしてみましょう。コードは以下の通りです。

let express = require("express")
let app = express()

app.get('/', (req, res)=>{
    res.send("Hello, World")
})

app.listen(process.env.PORT, ()=>{
    console.log("Server listening")
})

「Hello,World」と表示されたでしょうか。静的な文字ばかり表示していても面白くないので、今度はHTMLファイルを表示させましょう。まず、viewsフォルダにindex.htmlを作ります。ターミナルに

mkdir views
touch views/index.html

とコマンドします。そしてserver.jsに戻って

let express = require("express")
let app = express()

app.get('/', (req, res)=>{
    res.sendFile(__dirname + '/views/index.html') // この行を変更
})

app.listen(process.env.PORT, ()=>{
    console.log("Server listening")
})

とします。そしてviews/index.htmlに

<!DOCTYPE html>
<html>
    <head>
        
    </head>
    <body>
        <p><em>Hello. World!</em></p>
    </body>
</html>

とします。ではRunしてみましょう。斜体のハローワールドが表示されましたでしょうか。　ここまで見てきてわかる通り、鯖立てるのめっちゃ簡単ですよね。合計20行にも満たないコードですが、HTML表示までできました。ここから肉付けしていきましょう。いよいよSocket.ioの導入です。

　Socket.ioの導入

この後どんどんクライアントサイドのjsも書いていくのですが、いったんソケットの話をしましょう。わかりやすくするために図を作りました。

例えば、クライアントAがソケット通信でイベントを発生させると、それが鯖に通知され、それを処理します。

鯖は処理した内容に従ってまた新たなイベント（今回は同じイベントを発生させた）を発生させます。

鯖からイベント通知を受け取ったクライアントB、Cはイベントを処理します。

このような流れでリアルタイム通信が成り立ちます。以上のことを踏まえてソケット通信をしてみましょう。コードを書く前に、

mkdir js
touch js/index.js

とコマンドを打ってindex.jsを作っておきましょう。まずはサーバサイドから書いていきます。

let express = require('express')
let http = require('http')
let app = express()
let server = http.createServer(app)
let io = require('socket.io').listen(server);
  
app.use(express.static(__dirname))

server.listen(process.env.PORT);

app.get('/', (req, res)=> {
    res.sendFile(__dirname + '/views/index.html');
});

io.sockets.on('connection', (socket)=>{
    socket.on('eventA', (eventData)=>{
        io.sockets.emit('eventB', {msg: eventData.message})
    })
})

上のrequire系はもう呪文だと思って書いています。ここで

app.use(express.static(__dirname))

ですが、これを行うことでjsファイルやcssファイルを読み込むことができます。最初読み込んでるはずのjsファイルが見つからないと言われ私はSANチェックだったのですが、ドキュメントを見てみると、ミドルウェアという仕組みを使って静的ファイルを読み込むことで解決するとのことでした。そしてsocketの部分である

io.sockets.on('connection', (socket)=>{
    socket.on('eventA', (eventData)=>{
        io.sockets.emit('eventB', {msg: eventData.message})
    })
})

この部分ですが、簡単に言うと、eventAというイベントが発生したら、eventBというイベントを発生させよ、と言っているだけです。でもってeventAからはmessageというインデックスを持ったハッシュが渡されているわけです。このハッシュをeventDataという引数に格納して即時関数を実行し、新たにmsgというインデックスを持ったハッシュを渡しています。クライアントサイドも見てみましょう。HTMLのほうはjqueryとsocket.ioと先ほど作ったjs/index.jsを読み込み、ボタンにonclickイベントを付け加えただけです。

<!DOCTYPE html>
<html>
    <head>
        
    </head>
    <body>
        <p><em>Hello. World!</em></p>
        <input type="text" name="msg-input" id="msg-input"/>
        <button id="msg-btn" onclick="sendEvent()">Send</button>
        <div id="msg-whole"></div>
        
        <!--追加-->
        <script src="//ajax.googleapis.com/ajax/libs/jquery/1.8.0/jquery.min.js"></script>
        <script src="/socket.io/socket.io.js"></script>
        <script type="text/javascript" src="js/index.js"></script>
    </body>
</html>

var socketio = io()

var sendEvent = function(){
    msg = $('input#msg-input').val()
    socketio.emit('eventA', {message: msg})
}

socketio.on('eventB', function(data){
    $('div#msg-whole').prepend("<div>"+ data.msg +"</div>")
})

さて、クライアントjsですが、至極シンプルなコードですね。ボタンを押したらsendEvent()関数が呼び出されるのですが、そこではinputボックスの中身をmessageというインデックスのハッシュに入れて渡しています。先ほどやりましたね。そしてeventAが発生すると繋がっているすべてのクライアントPCにeventBが発生するので、そのイベントハンドラを記述しています。

socketio.on('eventB', function(data){
    $('div#msg-whole').prepend("<div>"+ data.msg +"</div>")
})

この部分ですね。単純にmsg-wholeのidを持つdivにちっちゃいdivをどんどん詰め込んでるだけです。このようなコードから、 1. メッセージを送信 2. 鯖が受け取り、受け取ったメッセージをすべてのクライアントに送信 3. クライアントが受け取り、divの中に追加という構図が浮かびますでしょうか。実際にタブを複製してやってみてください。リアルタイムで通信しているのがよくわかります。

私はこんな感じで作りました。

さぁ、ソケットは理解できたでしょうか？そもそも私が初心者なので拙い説明になってしまったと思うので、わからないところがあったら別途調べていただければと思います。しかし、ここまで理解できれば、あとはフロントを整えたり、イベントを増やしたりするだけでシンプルなチャットはできます。ぜひやってみてください。その一例として、私が実際に書いたコードを載せて、少し補足をさせていただきます。

.
├── js/
│  ├── index.js
│  └── def events.js
├── views/
│  └── index.html
├── style/
│  ├── index.sass
│  ├── _vars.sass
│  ├── _baseDesign.sass
│  ├── _msgStyle.sass
│  └── index.css
├── node_modules/
├── server.js
├── package.json
└── README.md

package.json

{
  "name": "express-socket-test-app",
  "version": "0.1.0",
  "description": "",
  "main": "server.js",
  "scripts": {
    "test": "echo \"Error: no test specified\" && exit 1",
    "start": "node server.js"
  },
  "author": "drumath2237",
  "license": "ISC",
  "devDependencies": {
    "express": "^4.16.2",
    "socket.io": "^2.0.4"
  },
  "dependencies": {
    "express": "^4.16.3",
    "socket.io": "^2.0.4"
  }
}

server.js

これはあまり変わりませんね。変更点としては、メッセージを送ったユーザーの名前がわかるようにしたことです。こうすることによって、自分のメッセージと他人のメッセージによってスタイルを変えることができますので。

let express = require('express')
let http = require('http')
let app = express()
let server = http.createServer(app)
let io = require('socket.io').listen(server);
  
app.use(express.static(__dirname))

server.listen(process.env.PORT);

app.get('/', (req, res)=> {
    res.sendFile(__dirname + '/views/index.html');
});

io.sockets.on('connection', (socket)=>{
    socket.on('message', (data)=>{
        io.sockets.emit("message",
        {
            msg: data.message,
            name: data.name
        })
    })
})

views/index.html

これもあまり変えていませんね。

<!DOCTYPE html>
<html>
    <head>
        <link rel="stylesheet" href="style/index.css" type="text/css" />
        <link href="https://fonts.googleapis.com/css?family=Itim" rel="stylesheet">    </head>
    <body>
        <h1>Hello,<input type="text" id="user-name-input" placeholder="put your name">! Send Messages!</h1>
        <input type="text" name="msg-input" id="msg-input"/>
        <button id="msg-btn" onclick="sendMessage()"> Send </button>
        <div id="msg-whole"></div>
        
        <script src="//ajax.googleapis.com/ajax/libs/jquery/1.8.0/jquery.min.js"></script>
        <script src="/socket.io/socket.io.js"></script>
        <script type="text/javascript" src="js/index.js"></script>
        <script type="text/javascript" src="js/def events.js"></script>

    </body>
</html>

js/index.js

makeMessage関数が少しややこしくなっていますが、スタイルを整えるためにやむなく･･･。

var socket = io()

socket.on('message', function(data){
    addMessage(data.msg, data.name)
})

var sendMessage = function(){
    var msg = $('input#msg-input').val()
    var name = $('input#user-name-input').val()
    socket.emit('message',
    {
        message: msg,
        name: name
    })
    $('input#msg-input').val('')
}

var addMessage = function(msg, username){
    var msgWhole = $('div#msg-whole')
    msgWhole.prepend(makeMessage(msg, username))
}

var makeMessage = function(msg, username){
    if(username===$('input#user-name-input').val())
        return "<div style='text-align: right'><div class='my-msg'><span class='user-name'>" + username + ":</span><span class='message'>" + msg + "</span></div></div></br>"
    else
        return "<div style='text-align: left'><div class='other-msg'><span class='user-name'>" + username + ":</span><span class='message'>" + msg + "</span></div></div></br>"
}

js/def events.js

ボタンを押したときのアニメーションやエンターキーでMessageを送信するようなイベントハンドラを記述しています。jQueryしか使ってない。

$(function(){
    $('button#msg-btn').on('mousedown', function(){
        var button = $(this)
        button.css({
            'top': '3px',
            'border-bottom': '0'
        })
    })
    $('button#msg-btn').on('mouseup', function(){
        var button = $(this)
        button.css({
            'top': '0',
            'border-bottom': '3px solid white'
        })
    })
    
    $('input#msg-input').keypress(function(e){
        if(e.which == 13){
            sendMessage()
        }
    })
})

style/index.sass

さてSassでスタイルを整えていますが、htmlに反映するためにはコンパイルしてcssファイルにする必要があります。毎度毎度やっているのはめんどくさいので私の場合はターミナルをもう一個起動して

sass index.sass:index.css --style expanded --watch

としておけばSassファイルを監視してくれるので保存するたびに自動でコンパイルしてくれます。

@import "baseDesign"
@import "vars"
@import "msgStyle"

html, body
    margin: 0
    padding: 0
    background: white
    text-align: center
    background: $bg-color
    color: white
    font-family: 'Itim', cursive
    height: 100%
    input
        @include base-design
        border: 0px
        border-bottom: 2px solid white
        &#user-name-input
            width: 200px
            text-align: center
    button#msg-btn
        @include base-design
        border-radius: 3px
        height: 40px
        box-shadow: 0px 3px 0px 0px white
        position: relative
        top: 0
    div#msg-whole
        height: 50%
        width: 80%
        border: 2px white dashed
        margin: 0 auto
        margin-top: 20px
        overflow: auto
        padding: 10px
        div.my-msg
            @include msg-style
            background: #96ed9e
            color: white
            &:after
                content: ""
                position: absolute
                top: 20px
                right: -15px
                height: 0
                width: 0
                border-top: 10px transparent solid
                border-right: 0px
                border-bottom: 2px transparent solid
                border-left: 15px white solid
            span.user-name
                font-weight: bolder
                font-size: 18px
                margin-right: 10px
            span.message
        div.other-msg
            @include msg-style
            background: $bg-color
            color: white
            &:after
                content: ""
                position: absolute
                top: 20px
                left: -15px
                height: 0
                width: 0
                border-top: 10px transparent solid
                border-right: 15px white solid
                border-bottom: 2px transparent solid
                border-left: 0
            span.user-name
                font-weight: bolder
                font-size: 18px
                margin-right: 10px
            span.message

style/_vars.sass

変数などを定義するパーシャルですが、背景色しか定義しなかったんで全然varsみたいに複数形にすることなかった。

$bg-color: #96d7ed

style/_msgStyle

メッセージのスタイルを記述したパーシャルです。

@mixin msg-style
    border: 2px white solid
    line-height: 40px
    border-radius: 15px
    min-width: 100px
    max-width: 50%
    display: inline-block
    margin: 10px
    padding-left: 10px
    padding-right: 10px
    font-size: 20px
    position: relative
    word-wrap: break-word

style/_baseDesign.sass

コントロールのスタイルを規定したパーシャルです。

@mixin base-design
    border: white 2px solid
    line-height: 30px
    background: rgba(0,0,0,0)
    font-size: 30px
    color: white
    font-family: 'Itim', cursive
    padding: 5px
    margin: 5px
    &:focus
        outline: 0

herokuにデプロイだ！

Gitを整える

まずgit周辺をいじりましょう。コミットするだけなんですけどね。

git init
git add .
git commit -m "first commit"

とすれば、大丈夫です。

herokuにあげる

いよいよherokuにあげちゃいましょう。まずherokuのアカウントを作ってください。無料でクレジットカードなしだと５つまでしか上げられないみたいで、さっき怒られました。 herokuのアカウントを作成出来たらターミナルでherokuにログインします。

heroku login

メアドとパスを入力します。前はパス入力中は何も表示されなくて怖かったのですが、今はアスタリスクが表示されて安心します。ログインで来たら続けて以下のように打ちます。

heroku keys:add
heroku create
git push heroku master
heroku rename "<好きな名前（既に使われているのは使用不可）>"

pushは時間かかるので少し待ちます。さてではherokuのマイページへ移動します。renameでつけた名前をクリック→右上の「Open App」をクリックで表示できたでしょうか。お疲れ様です！無事デプロイできました！

おわりに

お疲れさまでした。初心者なりに記事を書いてみましたが、いかがだったでしょうか。私はサーバサイドに全然詳しくなくて、もちろん鯖も持っていなければNode.jsやRailsもthe素人ですので、この機会にちゃんと勉強したいなぁと思いました。それにしてもheroku便利ですよね。まさにブラウザだけで完結するWebアプリ開発。Cloud9もすごい便利なので（ソーシャルコーディングとかできるし）是非、この機会に使ってみてはいかがでしょうか。以上となります。最後まで読んでいただき、ありがとうございました！ Happy Coding!

2018-03-06

【小ネタ】区分求積法でπの近似値を求めるアルゴリズム

C++ 小ネタ

　どうも、drumathです。今回はちょっとしたチップスとして、πの近似値を求めたいと思います。使用する言語はC++です。

数Ⅲ生の十八番
区分求積法で近似する
コードを書く

数Ⅲ生の十八番

数Ⅲで積分をしていると、置換積分やらなんやらのおかげで一見関係なさそうなところで $\pi$ が出てきます。
今回はそんな積分の中から、受験生の十八番ともいえるやつを使います。

$\displaystyle \begin{align} \int_{0}^{1} \frac{1}{x^{2}+1} dx \end{align}$

について、 $x=\tan \theta$ と置くと

$\displaystyle \begin{align} x \to 0のとき、\theta \to 0、\\\ x \to 1のとき、\theta \to \frac{\pi}{4} \end{align}$

また、

$\displaystyle \begin{align} dx = \frac{1}{\cos^2 \theta}d\theta \end{align}$

となるので、

$\displaystyle \begin{align} \int_{0}^{1} \frac {1} {x^{2}+1} dx &=\int_{0}^{\frac {\pi} {4} } \frac{1}{\tan^{2} \theta + 1} \cdot \frac{1}{ \cos^{2} \theta} d\theta\\\ &= \int_{0}^{\frac {\pi} {4} } \frac{\cos^2 \theta}{\cos^2 \theta}d\theta\\\ &= \int_{0}^{ \frac {\pi} {4} } d \theta \\\ &= [ \theta ]_{0}^{\frac{\pi}{4}}\\\ &= \frac{\pi}{4} \end{align}$

となります。つまり、

$\displaystyle \begin{align} \pi = 4\int_{0}^{1} \frac{1}{x^2+1}dx \end{align}$

ということがわかりました。

区分求積法で近似する

ここで、区分求積法を用いると、

$\displaystyle \begin{align} 4\int_{0}^{1} \frac{1}{x^2+1}dx=4\lim_{n \to \infty} \frac{1}{n} \sum_{k=1}^{n} \frac{1}{(\frac{k}{n})^2+1} \end{align}$

となることから、 $n$ をだんだん大きくしていけば、 $\pi$ の近似値が求まるみたいです。

コードを書く

では先ほどの結果をもとにコードを書いてみました。

#include <stdio.h>
#include <math.h>

double pi_calc(double);

int main(void)
{
  for (int i = 1; i < 1000000; i++)
    printf("\r%d: %f", i, pi_calc((double)i));
  return 0;
}
// 区分求積法を用いたアルゴリズム
double pi_calc(double n)
{
  double res = 0;
  for (int k = 1; k <= n; k++)
    res += 1.0 / (1.0 + pow(((double)k / n), 2.0));
  return 4.0 * res / n;
}

$n$ を100万まで動かします。
コードでところどころdoubleを使っているのは、floatを使ってしまうと誤差補正のせいで値がぶれてしまうことを防ぐためです。

今回は以上になります。最後まで読んでいただき、ありがとうございました。

2018-02-28

レンダラー開発のための測光学覚書Vol.2　レイトレースのための基本的な物理量からレンダリング方程式まで

レン覚物理レンダリング

　どうも、前回に引き続き、『レンダラー開発のための測光学覚書』シリーズ第2弾です。前回は立体角についての数学の話だったので、今回からは物理の話に入っていきましょう。

基本的な物理量
マテリアルの決定
レンダリング方程式
最後に
- 参考文献

基本的な物理量

放射束

　光とは粒子であり、波である。これが今の物理の定説です。しかし未だに謎が多く、光とは何かというのは完璧に答えることはできません。しかし、光は私たちに世界を知覚させることができる、とても身近な存在ですよね。ゆえに経験からいろんなことがわかります。例えば、虫眼鏡で黒い紙に光を1点に集中させれば火が出るように、光はエネルギーを持っています。ここで紹介する「放射束」は単位時間あたりの光のエネルギーという定義で、測光の分野では最も基本的な物理量です。

$\displaystyle \begin{align} \Phi = \frac{dQ}{dt} \end{align}$

ここで、 $Q$ は光のエネルギー、 $t$ は時間です。

放射照度

単位面積当たりの放射束を放射照度といいます。

$\displaystyle \begin{align} E(x)=\frac{d\Phi}{dA} \end{align}$

放射強度

ここで、前回の内容である、立体角の概念を導入します。ラジアンの立体版でしたね。忘れてしまった方はこちらをどうぞ。
drumath.hatenablog.com
放射強度は単位立体角当たりの放射束として与えられる物理量です。

$\displaystyle \begin{align} I(\vec{\omega})=\frac{d\Phi}{d\vec{\omega}} \end{align}$

ですので、逆に半球について積分してやれば、放射束が求められます。

$\displaystyle \begin{align} \Phi = \int_\Omega I(\vec{\omega}) d\vec{\omega} \end{align}$

放射輝度

　放射輝度は、単位投影面積、単位立体角あたりの放射束として定義されます。先に式を出すと、

$\displaystyle \begin{align} L(x, \vec{\omega})=\frac{d^2\Phi}{\cos \theta dAd\vec{\omega}} \end{align}$

$\cos \theta$ がいきなり出てきましたが、これは定義が単位面積ではなく、「単位投影面積」ということに起因します。
例えば、下の図のように、光が地面に向かって $\theta$ の角度で入射しているとしましょう。
f:id:drumath:20180228002606p:plain

このとき、光が当たっている面積を $A$ とします。このときの放射照度は $\displaystyle E=\frac{\Phi}{A}$ ですね。
ここで、光の進行方向に垂直な面を考えます。すると、その面と光が当たっている面積は、図に示した通り、 $A\cos \theta$ となります。これを、投影面積といいます。
放射輝度は積分すれば、いろんな物理量が導けるので、シミュレーションでは大変重宝します。例えば放射束を出したい時には

$\displaystyle \begin{align} \Phi = \int_A \int_\Omega L(x,\vec{\omega})(\vec{n} \cdot \vec{\omega})d\vec{\omega}dx \end{align}$

上の式では、地面に対する法線ベクトルの単位ベクトルと、入射ベクトルの単位ベクトルの内積を取ることで入射角を求めることができます。この $(\vec{n} \cdot \vec{\omega})$ の項をコサイン項と呼びます。

マテリアルの決定

　ここまで光について見てきました。ここでは実際にレンダリングするとき、もう一つ重要な「どんな材質に光が当たったか」を決めるための関数をご紹介します。

BSSRDF（英：bidirectional scattering surface reflectance distribution function、日：双方向散乱面反射率分布関数）

　光が物質の表面に当たったら、入射角と同じ大きさの反射角で反射し、屈折の法則に従って屈折することは高校生、いや中学生で習いました。しかし現実はそんなに甘くありません。なぜなら完全に平らな面など存在しないからです。実際にはすごく小さな凸凹があって、その凸凹に当たって光はあらゆる方向に反射していきます。これを散乱というのでしたね。また、光はどの物質でも通過します。つまり物質の表面で屈折した光はそのうち物質中で反射して出てくるかもしれないし、反射を繰り返して物質を通過するかもしれないのです。このように物質を通過した光が反射して出ていくことを「表面下散乱」といいます。このように表面化散乱の様子を示す関数が、BSSRDFと呼ばれる関数です。

f:id:drumath:20180228153043p:plain — https://news.mynavi.jp/article/graphics-59/　から引用

上の画像からもわかる通り、人の肌の質感などのレンダリングではこの技術が用いられます。しかし、ちょっと複雑ですね。なので普通はBSSRDFを簡略化したBSDFを用います。

BRDF（英：bidirectional reflectance distribution function、日：双方向反射率分布関数）

　BRDFはBSSRDFを簡略化したモデルです。正確に言うと、BSSRDFの簡略版は、BRDFとBTDFを足しあさせたBSDFというモデルなのですが、拡散反射や鏡面反射だけのモデルの場合、BRDFで事足りるということなんだと思います。BRDFは入射位置からでる反射光をモデル化したもので、先ほどの表面下散乱を考えないもので、入射光による放射照度と反射光による放射輝度の比で与えられます。
ここで、式の中の物理量に以下のものを使います。（なんか今更な感じもしますが）

物理量	意味
$\Phi$	放射束
$E$	放射照度
$L$	放射輝度
$f_{r}$	BRDF
$\theta_{i}$	入射角
$\vec{\omega}$	入射光、及び反射光の方向
$\vec{n}$	入射面に対する法線

添え字の $i$ や $r$ はそれぞれ入射光関係(incidence)、反射光関係(reflectance)という意味です。

$\displaystyle \begin{align} BRDF: f_{r} (x, \vec { \omega }_{i}, \vec{\omega}_{r}) &= \frac{dL_{r}(x,\vec{\omega}_{r})}{dE_{i}(x, \vec{\omega}_{i})} \\\ &= \frac{dL_{r}(x, \vec{\omega}_{r})}{L_{i}(x, \vec{\omega}_{i})(\vec{n} \cdot \vec{\omega}_{i})d\vec{\omega}_{i}} \end{align}$

さて式変形についての説明ですが、

$\displaystyle \begin{align} L_{i}(x, \vec{\omega}_{i}) = \frac{d^{2}\Phi}{\cos \theta_{i} dAd\vec{\omega}_{i}} = \frac{dE(x, \vec{\omega}_{i})}{(\vec{n} \cdot \vec{\omega}_{i})d\vec{\omega}_{i}} \end{align}$

となるためですね。

BRDFの意味（高校生なりの解釈）

　私は最初BRDFの式を見て、なんでこうなるんだかわかりませんでしたw　というのも、なぜ放射照度と放射輝度という二つの違う量の比で表すのかが意味不明でした。しかもこんなに複雑な式で、さも微積させる気しかないような雰囲気ですよね。頭の中ではグラサンかけてる厳ついおっさんが突然出てきたみたいなイメージでした。絡みづらいわ。
しかし、BRDFもこんな風に変形すれば、目的は一目瞭然です。

$\displaystyle \begin{align} dL_{r} ( x, \vec { \omega }_{r} ) &= f_{r} ( x, \vec { \omega }_{i}, \vec{\omega}_{r} ) L_{i} (x, \vec { \omega }_{i}) (\vec{n} \cdot \vec { \omega }_{i}) d \vec { \omega }_{i} \\\ \Leftrightarrow L_{r} ( x, \vec { \omega }_{r} ) &= \int_{\Omega} f_{r} ( x, \vec { \omega }_{i}, \vec{\omega}_{r} ) L_{i} (x, \vec { \omega }_{i}) (\vec{n} \cdot \vec { \omega }_{i}) d \vec { \omega }_{i} \tag{1} \end{align}$

つまり何が言いたいかというと、BRDFなんてものは関数というよりもよりも、係数みたいなもんじゃないのか。
このようにすると、後述するレンダリング方程式とほぼ形は一緒です。つまり、いま私たちが求めたいのは放射輝度である。前述のとおり放射輝度さえわかればほかの物理量は積分すれば求まりますからね。だから普通はBRDFを求めよう、なんてことはしないんじゃないかな。
もう一回式(1)を見てみましょう。

$\displaystyle \begin{align} L_{r} ( x, \vec { \omega }_{r} ) &= \int_{\Omega} f_{r} ( x, \vec { \omega }_{i}, \vec{\omega}_{r} ) L_{i} (x, \vec { \omega }_{i}) (\vec{n} \cdot \vec { \omega }_{i}) d \vec { \omega }_{i} \tag{1} \end{align}$

目的は $L_{r}$ を求めることでした。被積分関数を見てみると、3つの項があります。左から順に、BRDF、入射光の放射輝度、コサイン項です。反射光の放射輝度を求めるのには、入射光の放射輝度が必要だということは、納得がいくと思います。またコサイン項があることによって放射輝度が一緒でも、入射角によって反射光の放射輝度が違うことも読み取れます。だから、BRDFというのは、入射光がどれだけ変化するかを表しているだけなのです。ここまでよろしいでしょうか。
半球で積分をしている意味は、面に対して、あらゆる方向から入射してくる光の和という意味ですね。このおかげで、グローバルイルミネーションが実現できます。

レンダリング方程式

さて、ついにレンダリング方程式です。とはいえ、先ほどの内容がわかっていれば余裕です。ようは出ていく放射輝度を求める方程式です。
出射する光というのは、表面から放射される光と反射される光の和で求まります。これは感覚的ですね。つまり、

$L_{r}$	反射光(reflected)の放射輝度
$L_{e}$	放射光(emitted)の放射輝度
$L_{o}$	出射光(outgoing)の放射輝度

とすれば、求める $L_{o}$ は

$\displaystyle \begin{align} L_{o} = L_{e} + L_{r} \end{align}$

となります。また、式(1)より

$\displaystyle \begin{align} L_{o}(x, \vec{\omega}_{r}) = L_{e}(x, \vec{\omega}_{r}) + \int_{\Omega} f_{r} ( x, \vec { \omega }_{i}, \vec{\omega}_{r} ) L_{i} (x, \vec { \omega }_{i}) (\vec{n} \cdot \vec { \omega }_{i}) d \vec { \omega }_{i} \end{align}$

と書き換えられます。お疲れさまでした、レンダリング方程式の出来上がりです！

最後に

ここまで学習するのにはとても時間がかかりました。しかし、レンダリングの原理が少しわかったような気がして、自信が付きました。既存のレンダラーなど見ていると、いつもは何気なく見ている言葉も、「これ知ってる！」ってなってなんだか楽しいですw　BlenderのCyclesレンダラーも、マテリアルの名前は「～～BSDF」という風になってますけど、材質決定の関数だなんて思いませんでした。これからも頻度は高くないかもしれないけど、このシリーズは続けていこうと思います。最後まで見てくださり、ありがとうございました。

参考文献

qiita.com
↑いつもお世話になっています。より詳しい内容が知りたい方はこのリンクへ飛んでみてください。
raytracing.hatenablog.com
↑教育用のレンダラーを作ってくださっています。日本語のコメントがあるって素敵
qiita.com
↑こちらもわかりやすく解説してくださっています。

2018-02-20

Blenderで水の中の表現

Yafaray レンダリング 3DCG Blender 海洋

　どうも、今回はYafarayレンダラーを使って海の中を作ってみます。あまり長くはならないと思うのでお付き合いください。

海のオブジェクトとその他の環境を作る
海洋の設定
レンダリングしてみる
おわりに
追記（2018/3/2）

海のオブジェクトとその他の環境を作る

　では早速作っていきます。まず、画面に平面のオブジェクトを作成し、「海洋モディファイア」を付与します。 f:id:drumath:20180220180251p:plain
すると、こんな感じになります。
f:id:drumath:20180220180313p:plain
そして海洋を囲むように空間を作ります。
私は上面を除去したCubeオブジェクトに「厚み付けモディファイア」を付与しました。別に上面が閉じていてもいいのですが、Lampでサンを使いたいので私は開けました。
今の画面はこんな感じ。中にカメラも入れてください。
f:id:drumath:20180220180633p:plain f:id:drumath:20180220180640p:plain

海洋の設定

　先ほど追加した「海洋もモディファイア」ですが、私の設定はこんな感じ
f:id:drumath:20180220180911p:plain
時間を変えると波が変化するので、いい感じの時間を選んでください。イイ感じの時間とは、波がまぁまぁ立っているときですかね。波がないと、海洋独特のコースティクスが生まれませんので。
　マテリアルも編集していきます。
f:id:drumath:20180220181212p:plain
IOR（屈折率）は水の値を設定しました。

レンダリングしてみる

　ではフォトンマッピングでレンダリングします。主要な設定は以下の通り。

名前	値
Lighting Method	Photon Mapping
深度	16
バウンス	9
サンプル	16

f:id:drumath:20180220182059p:plain

海洋モディファイアの解像度の値を7→10に変えてみましょう。
f:id:drumath:20180220182527p:plain

おわりに

今回はすごく早く書き終わりましたw　とはいえCyclesでの再現で困っていたのでYafarayでできてよかったです。これからもYafarayの記事を書くかもしれないのでよろしくお願いします。
Yafarayについての他の記事↓
drumath.hatenablog.com

追記（2018/3/2）

動画で作ってみました。コンポジットし忘れたので、AviUtlでグロー効果をかけました。
youtu.be

2018-02-18

レンダラー開発のための測光学覚書Vol.１　立体角について

レン覚 3DCG レンダリング物理

　どうも。大学生の勉強の予習になればと、レンダラー制作に向けて少しずつ勉強を始めている、drumathです。実際わからないことだらけなので、『レンダラー開発のための測光学覚書』というシリーズで、学習してきたことを整理していきたいと思います。もしよろしければ、解釈や式の導出などが間違っていた場合に教えていただけるとありがたいです。
　今回は本格的に光学の分野には触れませんが、物理量の基本となる立体角についてまとめていきます。

立体角とは？
- ラジアンの定義を振り返る
- 平面から立体へ
微小立体角を求める
- 設定
- 微分だ！
終わりに
- 参考文献

立体角とは？

ラジアンの定義を振り返る

　高校になっていきなり導入された弧度法ですが、なぜ $\displaystyle 60^\circ = \frac{\pi}{3}$ となるのか、説明できますか？「簡単だよ！」という方は飛ばしていただいて構いません。
弧度法の定義は $\displaystyle \theta = \frac{l}{r}$ で、 $\displaystyle \theta$ は弧度法での角度、単位はラジアンです。 $\displaystyle l$ は $\displaystyle \theta$ での弧の長さ。 $\displaystyle r$ は半径です。つまり、円周の長さと半径の比であらわされる値なのです。とても直観的だと思いませんか？対して度数法の場合、一周を360個に分けた内、どれくらいを占めるかという、とてもあいまいな値です。一周が360°なんて一体だれが考えたのでしょうか。古代文明の60進法からですよね、たぶん。
　数学Ⅲを習うと、弧度法を使うメリットがわかります。詳しくはほかの記事に書いてあることを参照してほしいのですが、 $\displaystyle (\sin x)' = \cos x$ のような三角関数の微分は弧度法を使うことによって、きれいな形で表されます。個人的な解釈なのですが、三角関数自体、直角三角形の辺の比で表される関数なので、弧と半径の比という定義はなじみやすいのではないかと考えます。

平面から立体へ

　少し話は外れましたが、要は円周の長さと半径の比がラジアンなわけです。平面角では長さを使ったので、立体角では面積を使ってはいかかでしょうか。ということで立体角の定義は

${ \displaystyle \omega = \frac{A}{r^2} }$

となります。 $A$ は面積で、 $r$ は円の半径です。立体角の単位は「ステラジアン」といいます。
定義に沿うと、全球の立体角は

$\begin{align} \displaystyle \omega &= \frac{A}{r^2} \\\ &=\frac{4 \pi r^2}{r^2} \\\ &= 4 \pi \end{align}$

となります。下のようなイメージです。
f:id:drumath:20180218152709p:plain

微小立体角を求める

設定

　光の物理量は微分と積分のオンパレードです。ですので、ここで微小立体角、 $d \omega$ を決めたいと思います。まず先ほどの画像の領域について、横が $a$ 、縦が $b$ の微小長方形で、立体角が $\omega$ と定めます。
f:id:drumath:20180218153603p:plain
当然、

${ \displaystyle \omega = \frac{ab}{r^2} \tag{1} }$

となりますね。

微分だ！

さて、では $a$ 、 $b$ が微小な量だとします。そうすると立体角 $\omega$ も微小な量になり、長方形領域は点のように小さくなります。ここで領域の方向を以下のように定めます。
f:id:drumath:20180218155923p:plain
ではまず微小 $a$ について考えましょう。長さ $a$ である横の辺を含む円周の長さは弧度法の定義から求めることができるでしょう。弧度法での円周の長さ $l$ は

${ \displaystyle l = r\theta }$

と表わせます。また、今回求めようとしている円の半径は $r\sin \theta$ ですので、円周の長さは $2\pi r \sin \theta$ です。いったん図にして整理しましょう。
f:id:drumath:20180218220530p:plain

では、なぜ円周の長さは $2\pi r \sin \theta$ なのでしょうか。それは中心角が $\phi$ である扇形の弧の長さが $r\phi \sin \theta$ と表わされ、今は $\phi = 2\pi$ の場合を考えているからです。そして、微小 $a$ は、この円周が微小 $\phi$ の場合を考えていることがわかります。よって

$\begin{align} \displaystyle a &= dl \\\ &= r\sin \theta d \phi \tag{2} \end{align}$

ということになります。
すみません、たぶんもっとうまい説明があるし、なにしろすごく遠回りな考えなのはわかっていますが、私が思考した結果がこれなのです。引き続き $b$ についても同じような導出を行います。このやり方が嫌だという人は読み飛ばしてください。申し訳ございません。
さて、続きを見ていきましょう。長さ $b$ である縦の辺を含む円は、半径 $r$ ですので、 $2\pi r$ となりますね。図にするとこんな感じ。
f:id:drumath:20180218222630p:plain

$a$ の長さを求めたときのように、この $2\pi r$ というのは $r\theta$ において $\theta = 2\pi$ だからでした。よって

$\begin{align} \displaystyle b &= dl \\\ &= rd \theta \tag{3} \end{align}$

やっと微小領域の面積が求まりますね。式(2)と式(3)を式(1)に代入します。すると、

$\begin{align} \displaystyle d\omega &= \frac{ab}{r^2} \\\ &= \frac{r \sin \theta d \phi \times r d \theta}{r^2} \\\ &= \sin \theta d \theta d \phi \end{align}$

となり、無事、微小立体角を導くことができました！

終わりに

今回は以上になります。立体角の概念は光学において基本的なものになりますが、私は学習当初あまりイメージできなかったので、数学Ⅲを習っている高校生にわかるように記事を書いたつもりです。最後まで読んでいただき。ありがとうございました。

参考文献

qiita.com

2018-02-06

OpenGLでobjファイルの3Dデータを表示してみる

3DCG C++ OpenGL Ruby

　こんにちは。今回は入門したてのOpenGLでBlender出力のobjファイルビューワを作りたいと思います！

作るものについて
- 実装プラン
- 開発環境など
いざ、Rubyサイドを実装だ
いざ、OpenGLで描画
- ファイル階層
- ObjTest.cppを書いていく
最後に

作るものについて

　今回は下の画像のような3Dビューワを作ります。↓

f:id:drumath:20180205232453p:plain — 自作モデルを読み込んだ3Dビューワ

ほんとにシンプルなビューワですね(笑)　見て分かる通り、テクスチャやマテリアルなどは複雑化しそうなので今回は実装しません。今回使用するのは頂点、辺、面のデータのみです。

実装プラン

　OpenGLをメインのライブラリとして使っていきますが、objファイル解析をRubyで実装しようと思います。その理由は、単純にRubyのほうがC++よりも文字列の処理がやりやすいからです。本当はオールC++でも全然できると思いますが、私が純粋なルビイストということもあり、Rubyを使用します。つまり、

objファイルをRubyで解析
無駄なデータを省き、簡単な形式にしてファイルを吐き出す
C++でファイル読み込み
OpenGLで描画

という流れで実行していきます。また、objファイルはBlenderによって出力された3Dモデルデータを使用します。

開発環境など

　使用する環境はこちらです。

項目	内容
OS	Windows10
ライブラリなど	OpenGL、freeglut
コンパイラ	Ruby2.2.4p230、bcc32
エディタ	Visual Studio Code 1.19.3
Blenderバージョン	Blender2.79

いざ、Rubyサイドを実装だ

Blenderでobjファイルにエクスポート

　例えば、BlenderでCubeを追加したときのデフォルトのオブジェクトをobjファイル出力してみましょう。

f:id:drumath:20180205235638p:plain — 一辺2の立方体

まず、[ファイル]→[エクスポート]→[Wavefront(.obj)]とクリックしていきます。そうすると、下の画面のようになります。
f:id:drumath:20180206000141p:plain
普通にエクスポートボタンを押したくなると思いますが、ここで「三角面化」を忘れないでください！今回は面が三角形であることを前提に頂点を結んでいきます。もし多角形の面があったらちょっと厄介なことになりますので、必ず設定してください。

objファイルの形式を知る

　先ほどの要領でobjファイルをエクスポートすると、このような内容が吐き出されます。

# Blender v2.79 (sub 0) OBJ File: ''
# www.blender.org
mtllib cube_obj_test.mtl
o Cube_Cube.001
v -1.000000 -1.000000 1.000000
v -1.000000 1.000000 1.000000
v -1.000000 -1.000000 -1.000000
v -1.000000 1.000000 -1.000000
v 1.000000 -1.000000 1.000000
v 1.000000 1.000000 1.000000
v 1.000000 -1.000000 -1.000000
v 1.000000 1.000000 -1.000000
vn -1.0000 0.0000 0.0000
vn 0.0000 0.0000 -1.0000
vn 1.0000 0.0000 0.0000
vn 0.0000 0.0000 1.0000
vn 0.0000 -1.0000 0.0000
vn 0.0000 1.0000 0.0000
usemtl None
s off
f 2//1 3//1 1//1
f 4//2 7//2 3//2
f 8//3 5//3 7//3
f 6//4 1//4 5//4
f 7//5 1//5 3//5
f 4//6 6//6 8//6
f 2//1 4//1 3//1
f 4//2 8//2 7//2
f 8//3 6//3 5//3
f 6//4 2//4 1//4
f 7//5 5//5 1//5
f 4//6 2//6 6//6

これは何を示しているのでしょうか。
　まず、

v -1.000000 -1.000000 1.000000

のような「v ～」ではじまる行は、頂点の座標を記録しています。Vはおそらく頂点の英語Vertexを意味しています。「v」につづく小数が左からx,y,z座標を表わしています。
そして、

f 2//1 3//1 1//1

のように「f～」ではじまる行は面のデータを記録しています。「f」のあとに「2//1」のような数字をスラッシュで区切ったような記述がみられますが、それは「頂点/テクスチャ座標/法線ベクトル」というデータを表わしています。今回使うのはい一つ目の数字だけですので、あとは無視してください。この1つ目の数字は、どの頂点を使うかを示しています。つまり「2//1」だった場合、一番左の「２」という数が、「2番目の頂点を使うぞ」と言っています。何番目の頂点かは「v」から始まる行の上から何番目かで決まります。
面の情報はどの三つの頂点を使うかを表しています。つまり

f 2//1 3//1 1//1

ならば、この面は２番目の頂点と3番目の頂点と1番目の頂点で構成されていることを示します。

Rubyでobjファイルコンバータを作る

　Rubyでは、objファイルを読み込み、あとでC++側でscanfするだけで簡単にデータが格納できるようなファイルを新たに作ります。というわけでコードはこんな感じ。

# vars
vertex = []
face = []

# const regexp
Face_regexp = /^f (\d+)\/\d*\/\d* (\d+)\/\d*\/\d* (\d+)\/\d*\/\d*/
Vertex_regexp = /^v (.*) (.*) (.*)/

# conversion obj file
File.open(ARGV[0], 'r') do |file|
  file.read.split("\n").each do |line|
    if line =~ Face_regexp
      face << [$1, $2, $3].map { |m| m.to_i - 1 }
    elsif line =~ Vertex_regexp
      vertex << [$1, $2, $3].map(&:to_f)
    end
  end
end

# output vertex data
File.open('vData.txt', 'w') do |file|
  vertex.each do |e|
    file.puts e.join(', ')
  end
end

# output face data
File.open('fData.txt', 'w') do |file|
  face.each do |e|
    file.puts e.join(', ')
  end
end

たとえば、さきほどの立方体のデータを読み込むと、

ruby obj_convert.rb cube.obj

vData.txtとfData.txtという二つのファイルが作成され、vData.txtには、

-1.0, -1.0, 1.0
-1.0, 1.0, 1.0
-1.0, -1.0, -1.0
-1.0, 1.0, -1.0
1.0, -1.0, 1.0
1.0, 1.0, 1.0
1.0, -1.0, -1.0
1.0, 1.0, -1.0

fData.txtには、

1, 2, 0
3, 6, 2
7, 4, 6
5, 0, 4
6, 0, 2
3, 5, 7
1, 3, 2
3, 7, 6
7, 5, 4
5, 1, 0
6, 4, 0
3, 1, 5

という内容が出力されていることがわかります。vDataのほうは、単純に座標のリストです。fDataのほうは、三角面を構成する座標のインデックスを表わしていますが、実際コードを見るとわかりますが、インデックスを-1しています。それは、objファイルは上から1,2,3...と数えますが、C++の配列では0,1,2...と数えるので、C++の仕様に合わせるためです。

いざ、OpenGLで描画

ファイル階層

　ファイル階層はこんな感じになります。

\GL
- freeglut.h
- freeglut_ext.h
- freeglut_std.h
- glut.h
- Makefile.am
- Makefile.in
\lib
- \x86
  - freeglut.dll
  - freeglut.exp
  - freeglut.iobj
  - freeglut.lib
  - freeglut.pdb
ObjTest.cpp
obj_trans.rb

ObjTest.cppを書いていく

　最初にコードを出します。私はこんな感じで書きました。

#include <windows.h>
#include <stdio.h>
#include "GL\freeglut.h"

#pragma comment (lib, "lib\x86\freeglut.lib")

#define ARRAY_MAX 10000000

float vertex[ARRAY_MAX];
int lines[ARRAY_MAX];
int vertexDataSize=0, lineDataSize=0;

void disp()
{
  glClearColor(0,0,0,1);
  glClear(GL_COLOR_BUFFER_BIT);

  glEnableClientState(GL_VERTEX_ARRAY);
  glVertexPointer(3,GL_FLOAT,0,vertex);

  // 頂点の描画
  glPointSize(3);
  glBegin(GL_POINTS);
  {
    for(int i=0;i<vertexDataSize;i++)
      glArrayElement(i);
  }
  glEnd();

  glColor3f(0.6,0.35,0);
  glBegin(GL_TRIANGLES);
  {
    for(int i=0;i<lineDataSize;i++){
      glArrayElement(lines[i*3]);
      glArrayElement(lines[i*3+1]);
      glArrayElement(lines[i*3+2]);
    }
  }
  glEnd();

  glColor3f(1.0,0.75,0);
  glBegin(GL_LINES);
  {
    for(int i=0;i<lineDataSize;i++){
      glArrayElement(lines[i*3]);
      glArrayElement(lines[i*3+1]);

      glArrayElement(lines[i*3+1]);      
      glArrayElement(lines[i*3+2]);

      glArrayElement(lines[i*3+2]);
      glArrayElement(lines[i*3]);
    }
  }
  glEnd();
  glFlush();
}

void timer(int x)
{
  glRotatef(1,0,1.0,0);
  glutPostRedisplay();
  glutTimerFunc(50,timer,0);
}

void InitialProc()
{
  FILE *fpVData, *fpFData;
  fpVData = fopen("vData.txt","r");
  fpFData = fopen("fData.txt", "r");

  if((fpVData==NULL)||(fpFData==NULL)){
    printf("file error!!\n");
    return;
  }

  while (fscanf(fpVData,"%f, %f, %f",&vertex[vertexDataSize*3],&vertex[vertexDataSize*3+1],&vertex[vertexDataSize*3+2])!=EOF)
    vertexDataSize+=1;

  while(fscanf(fpFData,"%d, %d, %d",&lines[lineDataSize*3],&lines[lineDataSize*3+1],&lines[lineDataSize*3+2])!=EOF)
    lineDataSize+=1;
}

int main(int argc, char **argv)
{
  InitialProc();

  glutInit(&argc, argv);
  glutInitWindowPosition(50,100);
  glutInitWindowSize(500,500);
  glutInitDisplayMode(GLUT_SINGLE|GLUT_RGBA);

  glutCreateWindow("Obj Test");

  glRotatef(30,1.0,0,0);
  glRotatef(30,0,1.0,0);
  glScalef(0.8,0.8,0.8);

  glutDisplayFunc(disp);
  glutTimerFunc(100,timer,0);

  glutMainLoop();

  return 0;
}

InitialProc()関数内ではvData.txtとfData.txtからデータを入力して配列に格納しています。こうすることで、あとでglArrayElementsで簡単に頂点を描画したり、辺を作ったりすることができます。
また、timer()関数で0.05秒ごとに回転させています。

最後に

今回はOpenGLの覚書ということで記事を書きました。まだまだ拙いコードではありますが、objファイルビューワが無事開発できてよかったです。最終的にはBlenderのような3DCGアプリをいつか作ってみたいです。