厳密な光線追跡 - その23

前回まででMathemaitcaパッケージの形に光線追跡の基本的な部分はとりあえずはできたんだけど、イマイチ気に入らない。いちおうオブジェクト指向風、というかMac OS XのCore Foundation風に書くことはできて、比較的シンプルなコードにはなった。

どこが気に入らないかと言うと、

1. まだシンプルさが足りなくて、ごちゃごちゃしているところがある
2. 実行効率が悪い

のふたつである。

シンプルさが足りないのは、いちばんは面形状オブジェクトの中身。コードを見てもらわないとわからないけど（コードはこの記事の真ん中へんにリンクがある）、面形状（コードの中ではshapeと呼んでいる）を平面、球面、2次曲面、非球面、その他、に分類してそれを指定してオブジェクト生成するようになっているが、その分類そのものはオブジェクトではなく、単にタグになっている。

これもオブジェクトとすべきだった。そのためには「継承」の概念があったほうがいい。まず一般的な形状というオブジェクトを作って、そのサブクラスとして具体的な形状を指定した方が簡単になる。

今回Core Foundationに似せて継承の機能をあえて持たせなかった。そのせいで記述が重複している感じが残ってしまった。

また、実行効率の点では、球面と平面のたった2面を追跡するだけで光線一本あたり1m秒ほどかかっていて、これでは実用には耐えない。

最初、オブジェクト指向風にすることの目的は、ひとつはコードをシンプルにしたい、ということと、内部状態を持って、繰り返し計算することはなるべくキャッシュするようにしよう、と思ったからだった。

コードを見てもらえばわかるんだけど、いちおう計算のための資源をオブジェクトが作られたときに内部に保持するようになっている。

しかし、実際には保持している式のほとんどはキャッシュの意味がなかった。それはMathematicaの評価メカニズムに従ったコードになっていない、というのが原因である。

つまり、計算のための内部作業は光線追跡の実行時にほとんどが行われるような保持の仕方になっている。

例えば光線追跡の最上位の関数はsequenceOfMediaオブジェクトが

sequenceMedia[traceRay] = 
  Function[{r}, 
   Apply[raySeries, ComposeList[sequenceMedia[traceList], r]]];

というような形で保持している。ところがこの無名関数では関数本体はこの時点では評価されず、ComposeList（Mathematicaの合成関数（f(g(h(x)))みたいなの）を作る関数で、一本の光線の追跡は各面での追跡する関数の合成関数をつくることに対応するため、こうなっている）の状態のまま保持されている。オブジェクトを生成する時点でこの合成関数は必要な数値がすべてそろっていて、数値関数に展開できる。そうすれば追跡を実行するときの効率はあがる。

しかし、そのためには数値関数に展開できるようにそれぞれの中身も揃ってないといけない。しかし全部をそのように書くことができなかった。

この問題は、特に面形状オブジェクトを設計し直せば解決できる。しかし、そこまでやるか、と言う気がしてしまった。詳しいことは割愛するが、どちらしにても、光線追跡のような特定問題の計算では、Mathematicaはプロトタイプのためと割り切った方がいい、と思い直した。

Mathematicaで書けば表現の自由度は高いし、数学の式の通りに書けるので実装の手間はかなり軽い。それでちゃんと動作するものを作っておいて、実際の計算にはObjective-CやC++などに展開して、Mathematicaはオブジェクト設計の正しさの確認や、計算の答え合わせに使う、というのが、単にコードの実行効率だけでなく、実装の効率も含めた全体の効率向上（実際に光線追跡をそのあとどのくらい実行するか、にもよるけど）につながる、と考えた。

5.3  ということで・・・

ということで、とりあえず今回のMathematicaコードは

ここ

に置いておく。ついでにHoyaとSchottとOharaのガラスカタログを読み込んで波長に対する屈折率を計算するオブジェクトを作るパッケージも含めた。

しかし、これを使って光線追跡をしよう、とは考えないで欲しい。これは参考、あるいは中途半端なMathematicaコードの「反面教師」だと思っていただきたい。

また、追跡のための式はまとまっているほうが読みやすいので、今回の一連の記事をLaTeXでpdfにしたものを同梱した。これが光線追跡のスタンダードの考え方ではないし、ブログの記事にするためにかなり冗長にはなっているけど、古いFORTRANではなく、C++やJavaなどで近代的な光線追跡を実装しようと思う人には便利なような数学をまとめたつもりである。ただし、幾何学が苦手な人のためには書いていないのでベクトルの基礎などは補っていただきたい。

そして、Mathematicaによる光線追跡エンジンは今回でおしまいにして、これをもとにObjective-Cで実装する、というのを始めることにする。

乞うご期待。

5.4  実行テスト

最後にどんなものができたか、を紹介しておく。

とりあえずほんとうに光線追跡エンジン本体のみ、つまり連続した媒質と面形状の組に対して光線を与えると、最終面までの光線を追跡する、というもの。

例えば、平凸レンズの追跡をしてみる。

下の図は2枚の面をK-SFLD1ガラス硝材ではさんだものを作ってみたところ。

createSurfaceObjectは面形状オブジェクトを生成して、createBoundaryObjectはそれを含んだ面の位置や向きを保持するバウンダリオブジェクトを生成している。球面ひとつと平面ふたつのみっつのバウンダリオブジェクトを生成して、createSequenceOfMediaObjectで光学系全体を生成する。第1面（球面）と第2面（平面）の間がガラスで、その周囲が真空（vacuum）である、と定義している。

これに対して光線を2次元的に入射させているのがこれである。

traceRayというのが引数の光線を追跡する関数である。追跡時間をTimingで測定していて、81本の光線追跡に約90m秒かかっている。

入力の16番は追跡結果のひとつを表示させている。raySeriesというコンテナがあって、それにrayオブジェクトが並んでいる。rayオブジェクトはひとつめの引数が出発点、ふたつめが進行方向を表していて、それぞれが均一な媒質内の直線の光線要素になっている。

追跡結果をグラフ表示させてみたのが次である。

光軸をzにとっているので、光線は下から上に走ることになる。ちなみに今回のエンジンにはまだ表示の機能を実装していないので、raySeriesToLineというユーティリティ関数を作って手動で描いている。

焦点付近で大きな球面収差が発生していることがわかる。幾何光学を知らない人にはなんのことだかわならないけど、まあ、近軸追跡では起こりえないいかにももっともらしい結果がでた、と思ってもらえばいい。

もう少し光線本数を増やしていわゆるスポットダイアグラムを描かせたのがこれである。

さっきと同じではつまらないので、入射光線をy軸周りに5°回転させている。大きな球面収差とすこしコマ収差がでているのがわかる。この表示もエンジンの機能ではなく、手動で描かせている。

光軸を含む断面内で追跡して縦収差を描いたのがこの図である。

とうぜん大きな負の球面収差が発生していることがわかる。しかしちょっとバグが残っていてメッセージがでている。

コメント 8

凄い…と感嘆しながら、コード追いかけ・勉強しています（コード量も凄いです）。Mathematica9で実行しましたが、glassCatalogs.nbのファイル読み込み部の冒頭に、
SetDirectory["hoge"]でカレントディレクトリ設定追加が必要だった以外、何の問題も無く動きました（Part::partd: "部分指定Null[[-1]]の長さはオブジェクトの深さを超えています"もちゃんといます）。
by jun hirabayashi (2013-02-24 21:46) 

コメントありがとうございます。
Mathematica9でもそのまま動きますでしょうか。
GlassCatalogs.nbの読み込み設定は、自分の環境でしかデバグしていなかったため忘れていました。データの入ったディレクトリが$Path上にあれば読み込めます。
今Objective-Cに移植したいので、いろいろ動かしてみていると、けっこう動作のおかしいところがいっぱい見つかりました。
ちょっとなさけないです。
こうすれば、もっとシンプルになる、というようなご提案がありましたら、ぜひお願いします。
by decafish (2013-02-25 08:08) 

Mathematica9でもそのまま動きます。それにしても、こういう「設計にあたっての捉え方・考え方」が詰まった話やコード、とてもありがたい限りです。
by jun hirabayashi (2013-02-26 05:43) 

そうおっしゃっていただけるとうれしいです。
誰の、何の役に立つのかわからない、と思っていたのですが、少しでもお役に立てたというのはよかったです。でももともとは、単なる未来の自分への備忘録だったので中身はいいかげんですが。
ところで、Mathematicaは5と6の間の非互換性が一番大きかったようですね。古いバージョンを使い続けないといけない者にとっては、とりあえずは精神衛生上よいことです。
by decafish (2013-02-26 07:43) 

コードサンプル、再公開していただけませんか。
最近, 数式処理システムで厳密な光線追跡ができるのでは？
と思いついてから、夜も眠れません。
私としては、sympyに移植してみたいです
by msfjs (2017-06-20 00:13) 

ファイルの置き場所にDropBoxを使っているのですが、仕様が変わって直接リンクができなくなったようです。記事中のリンクをDropBoxのダウンロードページへのリンクに修正しました。そこからダウンロードしてください。

でもちょっと言い訳をさせてもらいますと、これはどっちかというと僕がObjective-Cにハードコードするためのプロトタイプとして作ったもので、最終的には数値計算をしていて数式処理としてのメリットを発揮しているとは言えないですし、古くてMathematica6.0前提のコードのままなので、ご期待にそえるかどうかよくわかりません。

symPyを使ったことが僕はないのですが、任意精度数、無限精度数やグレブナー基底が実装されているようなので代数演算は問題ないのですね。こういうのをPythonで実装するのってすごいです。どうしてもLISP型のLISTを効率的に処理できないと実用的な実装は難しいのではないかと思い込んでいました。面白いです。

移植できたら是非様子をお教えください。Mathematicaより速いかもしれません。
よろしくお願いします。
by decafish (2017-06-20 18:16) 

リンクの修正ありがとうございます。感謝します。

厳密な光線追跡の計算方法を学ぶことで基礎から光学を理解したい、数式を手計算で追いかけるのは面倒だけど、(無料で使える)数式処理システムとしてsympyを見つけた、
というのが本音です。

codeV10.8でABCD行列が導入され、
off-Axial光学系の近軸量にあたるものが計算できるようになりましたが、
三次収差や高次収差に当たるものが計算できれば面白いかなと
思っております。

進展がありましたら連絡します。
他の記事もたいへん勉強になるので、読ませていただきたく。
今後ともよろしくお願いします

ありがとうございました。
by msfjs (2017-06-20 22:45) 

お役に立てるのなら光栄です。

off axisで近軸量が定義できるんですか？とりあえず一本の光線を基準にするということなんでしょうか。近軸で収差だらけにする人はいないでしょうから、意志を持ってここを光軸にする、と宣言するという感じでしょうか。

symPy今まで知らなかったのですが、すごく充実しているようですね。驚きました。Wikipediaに「13000行」とあったので、まだ基本的なところだけかと思ってたんですが、Wikipediaは一桁二桁間違っているとしか思えません。

面白そうな話がありましたらまたお教えください。
よろしくお願いします。
by decafish (2017-06-21 10:02)