NO_WAIT

主にプログラミング

Penroseの部屋の鏡面反射をJavaScriptでシミュレート

Illumination ProblemまたはUnilluminable Room Problemと呼ばれる問題があります。 鏡張りの部屋の中に1点の光源を置いたとき、光の届かない一角は生じ得るかという問題です。 または、そのような一角をもつ部屋をどう作るかという問題とも考えられます。 若きRoger Penroseが考案した奇妙な形の部屋が目に付き、 プログラム(JavaScript)によるシミュレーションをしてみたくなりました。 以下、それを実装したものです。

Penroseの部屋

Canvas領域上でクリックすると光源の位置を変更できます。 光源から光線がオレンジ色の線として発生します。これらは部屋の壁で鏡面反射します。 Show Foci、Show Light Sourceをチェックするとそれぞれ上下の大きな楕円の焦点を赤い点で、 光源を青い点で表示します。

このPenroseによる不思議な形の部屋はその中のどこに光源を置いても、 光の届かない場所が必ずできるという性質があります。

f:id:shinaisan:20140617223618p:plain

「きのこ」の裏が両側とも影になる場合

Penroseの部屋の奇妙な特徴である、両側面から内側に向かって生えたきのこ状の突起物に注目します。 きのこの傘の両端は大きな楕円の焦点に接しています。 ここで部屋の中腹部のあたりに光源を置いてみます。 具体的には、4つの焦点からなる長方形の内部にあり、左右のきのこ傘で挟まれた領域に光源を置きます。 すると図のように、傘の裏側にある4つの矩形には光が侵入しなくなります。

f:id:shinaisan:20140617223716p:plain

「きのこ」の裏のうち片側だけが影になる場合

上の例で置いた光源を、4焦点がなす長方形より上または下に移動してみますと、 きのこにより影になる領域は、光源の反対側のみになります。

f:id:shinaisan:20140617223739p:plain

f:id:shinaisan:20140617223751p:plain

部屋の上(下)半分のみを照らす

部屋の隅にある小さな矩形の一角に光源を置くとどうなるでしょうか。 きのこ傘の影にあたる部分です。 上側の影に置けば、部屋の上部のみが照らされ、上部楕円の焦点より下側には 光が達しなくなります。

f:id:shinaisan:20140617224324p:plain

楕円の部屋

上のプログラムで、Room TypeをElliptic Roomに変更すると、部屋が横長の楕円に変わります。 Single Rayチェックボックスで光線を1本だけにすると、 楕円の部屋での光の反射状況がわかりやすくなります。 光線は焦点を囲むような反射を繰り返すか、焦点の間(焦点を結ぶ線分上の点)を通る反射を繰り返すかのどちらかになります。

f:id:shinaisan:20140617224345p:plain

f:id:shinaisan:20140617224408p:plain

ちょうど光源を楕円の焦点の位置に置くと、すべての光線が 反対側の焦点を目指すようになります。 焦点から楕円上の任意の点で反射した光が別の焦点へ達するこの性質はよく知られているようです。

f:id:shinaisan:20140617224431p:plain

結局これら楕円内での反射の性質がPenroseの部屋の暗がりの存在を説明します。 例えば上で述べた「きのこ」の裏から発した光は焦点の周りを通るので、 焦点を結ぶ線分上の点を通らず、部屋の反対側へは光が届かなくなります。

楕円の性質を計算によって確かめる

焦点から発した光が別の焦点へ向かうことを数式処理システムSageを用いて確かめてみます。 (手計算で十分なのですが、単なる趣味です。 数式処理システムを少しは使えるようになりたかったので…)

長径\(2a\)、短径\(2b\)、原点中心の楕円を想定します。横長です。

\[ {x^2 \over a^2} + {y^2 \over b^2} = 1,\ a \geq b \gt 0 \]

上記楕円上の点\(P(x, y),\ y \ne 0\)における楕円の接線の傾きは正接で表すと

\[ t \stackrel{\rm def}{=} \tan{\theta} \stackrel{\rm def}{=} {dy\over dx} = - {{b^2 x} \over {a^2 y}} \]

となります。

焦点\((f, 0)\)から点\(P\)へ引いた直線の傾きも正接で以下のように表します。

\[ t_1 \stackrel{\rm def}{=} \tan{\theta_1} \stackrel{\rm def}{=} {y \over {x - f}} \]

もう一方の焦点\((-f, 0)\)についても同様にします。

\[ t_2 \stackrel{\rm def}{=} \tan{\theta_2} \stackrel{\rm def}{=} {y \over {x + f}} \]

目的は

\[ \tan{(\theta_1 - \theta)} = \tan{(\theta - \theta_2)} \]

を示すことです。

計算してみます。

var('a b f x y')
def tdiff(t, u): # 正接の加法公式: t = tan(x), u = tan(y)に対しtan(x - y)を計算する。
    return (t - u)/(1 + t*u)

t = - (b^2*x)/(a^2*y)
t1 = y/(x - f)
t2 = y/(x + f)

公式

\[ \tan{\left(\alpha - \beta\right)} = \frac{\tan\alpha - \tan\beta}{1 + {\tan\alpha}{\tan\beta}} \]

を関数tdiffとして定義しました。

\[ \tan{\left(\theta_1 - \theta\right)} - \tan{\left(\theta - \theta_2\right)} \]

を計算します。

tmp1 = tdiff(t1, t) - tdiff(t, t2)

楕円の方程式を用いて\(y^2\)を消去します。分子にしか興味がないのでnumeratorを取ります。 substituteではなくsubstitute_expressionを使います。

tmp2 = tmp1.numerator().substitute_expression(y^2 == b^2*(1 - x^2/a^2)).expand()

結果tmp2は次のようになります。

\[ -2 \, a^{4} b^{2} x + 2 \, a^{2} b^{4} x + 2 \, a^{2} b^{2} f^{2} x \]

さらに焦点の式\(f^2 = a^2 - b^2\)を用いて簡単にすると

tmp2.substitute_expression(f^2 == a^2 - b^2).expand()

0が得られ、求める結果が得られました。

楕円の性質の応用事例

上のような性質を応用した事例がないかと探していたら医療分野でありました。 ESWL - Extracorporeal ShockWave Lithotripsy (体外衝撃波結石破砕治療)という尿管結石の治療法がそれです。 衝撃波発生源位置と結石の位置を楕円の焦点に置くことで、楕円体内の反射衝撃波を結石に集中させ、これを破砕するという手法のようです。 (専門外なので間違っているかも…) 興味深い応用ではありますが、実際にお世話にはなりたくないですね。

Penroseの部屋 - 不完全版

ところで、プログラムにはRoom TypeとしてIncomplete Penrose Roomというものも用意しています。 これは「きのこ」の位置を焦点からランダムにずらしたもので、 この部屋においては光線がどこにでも侵入します。 Penroseの部屋において楕円の焦点が決定的な役割を果たしていることがわかります。

余談: Tokarskiの部屋

f:id:shinaisan:20140617225013p:plain

Room TypeにTokarski Roomという謎の部屋が用意されています。 実はPenroseの部屋の前にこのTokarski氏の考案したとされる多角形の部屋の存在を知ったのが プログラムを組んだ元々の動機でした。 深入りはしませんが、この部屋はある1点に光源を置くと、ある1点に光が到達しなくなるというものらしいです。 よくよく考えずともわかりますが、これをコンピュータで可視化して示すというのは不可能でしょう。 ここでいう1点とはどうやら1pixelのことでもなく、ある大きさの1区画というのでもない、「質点」であるからです。 これはvisualizeできないでしょう。 こんなことにも気がつかなかったとは何事でしょうか…

プログラムについて

プログラム中では線分と線分、線分と楕円の交点検出を行っていますが、 連立方程式2次方程式を愚直に解いているだけです。

困ってしまったのが、CanvasAPIに楕円の弧を描画する関数が見つけられなかったことです。 円の弧であればarcで描画できます。 仕方なく、次のようにscaleの変更で誤摩化しています。

    var h = this.radiusH / this.radiusV;
    ctx.save();
    ctx.translate(this.center.x, this.center.y);
    ctx.scale(h, 1);
    ctx.beginPath();
    ctx.strokeStyle = "black";
    ctx.lineWidth = 1;
    ctx.arc(0, 0, this.radiusV, this.startAngle, this.endAngle, false);
    ctx.stroke();
    ctx.restore();

もっと調べたいこと

  • より良い交点検出や反射の計算法。
  • 【募】Canvas APIで楕円の弧を描く方法
  • Tokarskiの部屋について。

参考