シェーディング言語で偏光顕微鏡のクロスニコル像の再現を試みた話

はじめに

こちらの記事は、エンジニア作業飲み集会 Advent Calendar 2022 10日目の記事です。

内容としては、2022/12/17, 18に開催されるバーチャル学会2022で発表する予定の内容に関する記事で、エンジニアとは直接関係ない内容になりますのでご了承ください。

また、内容的に間違いがないように努めていますが、ほぼ独学に近いので正確でない部分があるかもしれません。もし誤り等がありましたらTwitter（@cleantted_s）等でご指摘ください。

 

 

なんの話をするの？

ざっくりいうと、『「偏光顕微鏡」を使って見られる「クロスニコル像」というものを、shaderを使って再現できないか』というものです。

 

ひとつずつ解説していきます。

 

まず「偏光顕微鏡」ですが、簡単にいうと「観察する試料を挟み込むように2枚の偏光板を備えた顕微鏡」です。偏光した光を使って試料の光学的な性質を調べることができる顕微鏡です。地学での岩石や鉱物試料の観察の他にも、医療品や繊維の分野でも使われているそうです。

偏光顕微鏡は偏光板の状態によって大きく分けて2種類の像を観察することができ、「偏光板が1枚(もしくは2枚の偏光板が平行)の状態で観察する像」を「オープンニコル（開放ニコス）」、「偏光板が2枚で直交する状態で観察する像」を「クロスニコル（直交ニコル）」と呼びます。

オープンニコル像は肉眼で見るものとほとんど同じ見た目ですが、クロスニコル像では試料の光学的性質によって肉眼でみるものとは異なる見た目になります。どういう見た目になるかは、こちらのサイトやこちらの動画を見ていただければと思います（結構綺麗ですよ*1）。

 

www.ha.shotoku.ac.jp

www.youtube.com

 

ということで、このクロスニコル像をVRChat上での再現を、シェーダーを使ったらできそうだという着想から、シェーダー言語（GLSL）での実装を行ったというのが今回の話になります。

 

今回は「まず簡単な例での実現」を目指したため、いくつか制約*2を付けた状態に限っています。また、今回の内容は原理的には観察する対象によらず適用できるものなのですが、ユースケースとして「岩石・鉱物」での再現することを中心に話を進めていきます。

物理的な原理について

クロスニコル像を説明する上で重要になってくる現象は「複屈折」と呼ばれるものです。

複屈折で有名なのは方解石という鉱物です。方解石を通すことで反対にある文字が2重にみえ、方解石を回転させることで片方の像の位置が変化する様子が見られます。

youtu.be

ここでは詳しい複屈折の原理については触れませんが、ざっくりと説明すると光が通過するときに2つの速度の異なる光波に分かれるため、像が2つに分かれています。物体を通る速度が異なるということは、それぞれに対して異なる2つの屈折率を持つと解釈できます。

 

方解石を含め、一般的な結晶質の物質は方向によって性質が変わる「光学的異方体」となっています。反対に方向によって光学的な性質が変わらない物質を「光学的等方体」と呼びます。水やガラスなど非晶質な物質が多く見られます。「光学的異方体」に光が透過するとき、一般に互いに直交する振動方向に偏光した、2つの速度の異なる光波に分かれて進みます。

 

複屈折を示す物質であっても、方位によっては2つの光波の速度が一致する場合があります。このときの方位のことを「光軸(光学軸)」と呼んでいます。光学的異方体の結晶には光軸が1本の「一軸性結晶」と2本の「二軸性結晶」の2種類に大別できます。

 

ここで、「二軸性結晶のある方位についての、2つの光波に対する屈折率がどうなるか」ということを考えます。

二軸性結晶の屈折率は、ある互いに直交する3直線に関して対称になります。その直線（軸）を「主要振動軸（光学的弾性軸）」と呼び、X、Y、Zという記号で表されます。

3軸それぞれの方向での2つの屈折率はX軸で \beta と \gamma、Y軸で \alpha と \gamma、Z軸で \alpha と \beta になります（ここで \alpha、\beta、\gamma (\alpha \le \beta \le \gamma) は主屈折率と呼ばれる）。このとき、光軸はXZ平面に存在しており、2本の光軸はX軸またはZ軸に対して対称に存在します。ちなみに \alpha = \beta もしくは \beta = \gamma のとき、光軸は1本だけとなり一軸性結晶になります。

 

入射光の方位についての2つの屈折率 n_1、n_2 (n_1 \lt n_2) は、X、Y、Z軸とそれぞれ \alpha、\beta、\gamma で交わる三次元楕円体（「屈折率楕円体」と呼ばれています*3）を考えたとき、楕円体を入射光の方位に垂直な平面で切ったときにできる楕円の長軸と短軸の長さに一致します。入射光の方位が光軸と一致するとき、この断面の楕円は円になります*4。

二軸性結晶の屈折率楕円体（座標系は左手系）

ある入射光に対する屈折率 n_1、n_2 の解釈
オレンジ色の面は、入射光に垂直な面で切った屈折率楕円体の断面（楕円）

 

2つの屈折率 n_1、n_2 は以下の式で求められます。ここで、 \phi、\phi' はZの両側にある2つの光軸のそれぞれと、入射光の方位との間になす角度になります*5。

 

n_1^2 = \cfrac{2\gamma^2\alpha^2}{(\gamma^2+\alpha^2) + (\gamma^2 - \alpha^2)\cos{(\phi - \phi')}}

n_2^2 = \cfrac{2\gamma^2\alpha^2}{(\gamma^2+\alpha^2) + (\gamma^2 - \alpha^2)\cos{(\phi + \phi')}}

 

ところで、実は光軸の方位も主屈折率を用いて表すことができます。2つの光軸の間の角度の半分（つまり、1つの光軸とZ軸とのなす角度）を V_Z とすると、以下の関係が成り立ちます。

 

\sin{V_Z} = \cfrac{\gamma}{\beta}\sqrt{\cfrac{\beta^2 - \alpha^2}{\gamma^2 - \alpha^2}}

\cos{V_Z} = \cfrac{\alpha}{\beta}\sqrt{\cfrac{\gamma^2 - \beta^2}{\gamma^2 - \alpha^2}}

 

以上から、主屈折率 \alpha、\beta、\gamma を定めることによって光軸の方位が決定でき、屈折率を求めたい方位との角度を求めることができるため、任意の方位の屈折率 n_1、n_2 が求められます。

 

ここで、ようやく偏光顕微鏡で光学的異方体を観察することを考えます。

偏光顕微鏡では、光学的異方体の試料の上下に互いに直交する偏光板が2枚配置されています。前述した通り、光が試料を通過するときに2方向の光に分かれますが、入射光が偏光した光の場合、透過後の光の強度は（試料による光のエネルギーの吸収がないとすれば）入射光を各方向成分に分解したものになります。

オープンニコルの場合は、透過光をそのまま観察することになりますが、クロスニコルの場合、透過光が、もう一方の偏光板によって入射光と垂直な方向の成分だけが取り出されることになります。このとき、試料を透過する間に分かれた2つの光波には光路差があるため、波長によって干渉による強度に差ができます。この干渉によってできる像がクロスニコル像です。

 

入射光の振動方向と試料の透過後の光の振動方向のいずれかとが一致する場合、透過光は入射光と同じ振動方向だけのものになり、もう一方の偏光板で遮断されるため真っ暗になる。このような状況は、試料を入射光を軸に一回転させる間に4回あるため、一回転の間に明暗が4回繰り返すことになります。

 

具体的な透過光の強度を式で表します。

波長ごとの透過光の強度 I(\lambda) は、入射光の強度を I_0(\lambda) として、以下の式で求められます*6。

 

具体的な計算式方法は以下の通りです:

\bar{x}, \bar{y}, \bar{z} は CIE 1931 XYZ等色関数で、 Kは入射光での Y の値が 1/2.2 *7になるように定められた定数として、XYZ色空間の値を以下の式で求めます。

 

\begin{equation}
\left\{\begin{matrix}
 X = K \sum^{800}_{\lambda=400} \bar{x}(\lambda) I(\lambda) \\ 
 Y = K \sum^{800}_{\lambda=400} \bar{y}(\lambda) I(\lambda) \\
 Z = K \sum^{800}_{\lambda=400} \bar{z}(\lambda) I(\lambda)
\end{matrix}\right.
\end{equation}

 

さらに、XYZ色空間の値を以下の式でRGB色空間の値に変換します。もし、求められた値が0から1の範囲を超えていた場合は範囲に収まるように丸めています。

 

\begin{equation}
  \begin{pmatrix}
  R' \\
  G' \\
  B' \\
  \end{pmatrix}
  = 
  \begin{pmatrix}
  3.2410 & -1.5374 & -0.4986 \\
  -0.9692 & 1.8760 & 0.0416 \\
  0.0556 & -0.2040 & 1.0570 \\
  \end{pmatrix}
  \begin{pmatrix}
  X \\
  Y \\
  Z \\
  \end{pmatrix}
\end{equation}

 

\begin{equation}
  C =
  \begin{cases}
  12.92 C' & \text{ if } x \le 0.00304 \\
  1.055 C'^{1.0/2.4} - 0.055 & \text{ if } x > 0.00304 
  \end{cases}
  (C=R, G, B) 
\end{equation}

 

これで波長ごとの強度からRGBの値が求められるようになりました。

shaderでの実装

実際のshaderでの実装の話に入ります。

このあとの話で混乱を避けるため、 "グローバルな座標系" 、つまりUnity上で使われる座標系の軸は X'、 Y'、Z' と表記します。試料の主要振動軸のXYZ軸とは一般に一致しないので注意してください。主要振動軸の座標系は "ローカルな座標系" と読み替えるとわかりやすいかもしれないです*8。ベクトルの方位については特に断りがない場合、 "グローバルな座標系" で表現されるものとして解釈します。

 

さらに、話を簡単にするために、入射光の方位ベクトルは Y'軸正の方向と平行であるとして固定しておきます。このとき、入射光の振動方向は X'Z'平面に平行になります。

 

シェーダーに与える入力を整理すると、以下になります。試料を回転させたことによるクロスニコル像の変化も計算するため、試料の回転角度も入力に入れています。試料の回転軸は、入射光の方向ベクトルと同じ（つまり Y'軸）であるとします。

 

  • 試料の主屈折率: \alpha, \beta, \gamma
  • 試料の厚さ: d [nm]
  • 試料の初期姿勢（詳細は後述）
  • 試料の回転角度

試料の初期姿勢は、各ピクセルごとの X'、Y'、Z' 軸のオイラー角として RGB の値にエンコードしたテクスチャとして与えて、シェーダーでクォータニオンへと変換しています*9。

 

試料の初期姿勢を与えるテクスチャ

ここからクロスニコル像を計算するためには、以下のステップを踏む必要があります。

  1. 主屈折率 \alpha、\beta、\gamma からZ軸と光軸のなす角 V_Z を求める
  2. 試料の初期姿勢と回転角度と V_Z から、光軸の方位ベクトル（2つ）を求める
  3. 光軸の方位ベクトルと入射光の方位ベクトルが成す角 \phi、\phi' を求める
  4. 主屈折率と \phi、\phi' から、2つの屈折率 n_1、n_2 を求める
  5. 厚さ d と n_1、n_2 から、レターデーション R を求める
  6. 試料の初期姿勢から、初期姿勢での透過光の振動方向を求める
  7. 入射光の振動方向と透過光の振動方向とのなす角 \theta を求める
  8. 透過光の波長ごとの強度 I(\lambda) を求める
  9. 波長ごとの強度からRGBの値を求める

 

順番に必要な部分を補足していきます。

 

ステップ2の光軸の方位ベクトルを求める方法について、光軸が主要振動軸のXZ平面上にあることを利用すると、「Y軸を中心に、Z軸を V_Z （もしくは -V_Z）だけ回転させた方向」で求められることがわかります。ある軸を中心に回転させる操作はクォータニオンを用いることで簡単に表現できるため、初期姿勢をクォータニオンで表現しておくことで光軸の方位ベクトルが簡単に計算できます。

あとは、内積計算で2つのベクトルのなす角が求められます（ステップ3）。

 

ステップ6、7で \theta の値を求めるにあたって、「試料回転後の透過光の振動方向」を求める必要があります。ですが、透過光の振動方向は入射光と同じく X'Z'平面に平行であり、試料の回転軸は Y'軸であることから、「初期姿勢での透過光の振動方向」を求めておいて、後から回転させても同じ結果が得られることがわかります。なので、ステップ6では「初期姿勢での透過光の振動方向」を考えています。

「初期姿勢での透過光の振動方向」を求めるわけですが、執筆中に計算が間違っていることに気が付きました。なので、ここでは正しい求め方の気持ちだけを話しします*10。

「初期姿勢での透過光の振動方向」は、屈折率楕円体を光の入射方向に垂直な面で切ったときの断面にできる楕円の長軸・短軸と一致するということを前節で触れました。なので、初期姿勢から試料ではなく（相対的に）光軸の方位が変化していると解釈した上で、楕円体の切断面の楕円の長軸・短軸の向きを求めます。ここで、長軸・短軸の方位を r = (r_X, r_Y, r_Z)、s = (s_X, s_Y, s_Z) とすると、以下の式が成り立ちます。

 

\cfrac{r_X s_X}{\alpha^2} + \cfrac{r_Y s_Y}{\beta^2} + \cfrac{r_Z s_Z}{\gamma^2} = 0

 

一発でこの式が成り立つ r、s を見つけられればよいのですが、それはそれで複雑になるため、切断面上の適当な直交するベクトルを一旦 r、s として、光軸を中心にして回転させることで、上記式を満たすベクトルを探します。このときの回転角度 \omega は以下の式で求められます。

 

\omega = \cfrac{1}{2} arctan \cfrac{2 (\cfrac{r_X s_X}{\alpha^2} + \cfrac{r_Y s_Y}{\beta^2} + \cfrac{r_Z s_Z}{\gamma^2})}{\cfrac{(r_X^2 - s_X^2)}{\alpha^2} +\cfrac{(r_Y^2 - s_Y^2)}{\beta^2} + \cfrac{(r_Z^2 - s_Z^2)}{\gamma^2}}

 

これで透過光の振動方向が求まるので、入射光の振動方向とのなす角度を求めてあげて、試料の回転を反映させてあげれば \theta が求まります。

 

最後に、ステップ8、9ですが、実際の実装ではこの処理はシェーダーで行っていません。前節で説明した処理をシェーダーで実現しようとすると、非常に負荷がかかることが予想されます*11。なので計算時の負荷を少なくするように、以下のようにあらかじめ前計算を行って負荷軽減を図っています。

透過光の波長ごとの強度は、入射光の強度を固定すれば R と \theta の2変数のみで決まります。なので、事前に R を横軸に、 \theta *12 を縦軸に取ったときの干渉色の画像を作成しておき、これまでのステップで求めた R、\theta の値をもとにRGBの値を取得させています。

 

事前計算した干渉色の図
横軸の Rの最大値は 3000 [nm]として計算している

以上のステップを踏まえて、シェーダーを実装しました。実際の実装はGitHub上にリポジトリを上げているので、そちらを御覧ください。

github.com

最終的な青果物は動画を御覧ください。なお、オープンニコルが真っ白だとあまりおもしろくないのでテクスチャ追加しています。
最初の姿勢によって干渉色が変わっていることと、1/4回転の間に明暗を繰り返している様子がわかります。

youtu.be

最後に

かなり駆け足ですが、やったことを一通りまとめました。もし興味ある方はバーチャル学会当日にお越し下さい！

（VRChatのワールドURLは後日追記します）

*1:余談ですが、私は学生時代地質学が専門でしたが、研究していた石の偏光顕微鏡での見た目がきれいで好きだったので、研究の半分くらいのモチベーションが偏光顕微鏡で観察することでした。

*2:全体の屈折率は共通になる・オープンニコルでの鉱物色を考慮できていないなど

*3:「インディカトリックス」とも呼ぶようです。ちなみに、2つの屈折率がどう変わるかを示した「屈折率曲面」とは別物。

*4:都城・久城 (1972), 『岩石学Ⅰ』

*5:坪井 (1959), 『偏光顕微鏡』

*6:王子計測機器株式会社, 『直行ニコル干渉色とレターデーション』(URL)))。

 

I(\lambda) = I_0(\lambda)\sin^22\theta\sin^2{\cfrac{\pi R}{\lambda}}

 

ここで、 \theta は入射光の振動方向と透過光の速度が遅い方の光の振動方向（つまり n_2 の方向）との成す角度とします。 R は試料透過後の2つの光の光路差で、入射光に対する2つの屈折率を n_1、n_2 、試料の厚さを d として R = d(n_2 - n_1) で求められます。この R は「レターデーション」と呼ばれています。

 

以上から、無色の光学的異方体の場合、クロスニコル像は物体の主屈折率、物体の厚さと入射光に対する物体の方位によってのみ決まることがわかります。したがって、これらを入力として与えることでクロスニコル像を再現することが可能になります。

波長からRGBへの変換

前項で波長ごとの光の強度がわかったので、これを画像として落とし込むためにRGBの値に変換します。具体的には、「各波長の強度からXYZ色空間の値に変換し、その後RGB色空間の値に変換する」という方法をとります。日本語での解説はこちらの記事をご参照ください。

 

qiita.com

 

XYZ色空間の値をRGB色空間の値に変換する方法は、Stokesら(1996) をベースに行いました。CIE 1931 XYZ等色関数の値は「10-deg XYZ CMFs」の stepsize 1mmのものを用いています((値はここから取得: Colour matching functions

*7:正直、この値の妥当性の検証ができていませんが、それっぽい結果になるため一旦これで良しとしています。妥当性について誰かご指摘ください……。

*8:一般にUnity上は「左手系」ですが、主要振動軸の座標系は「右手系」で記載されることが多いです。Unity上で扱う上で主要振動軸も一応「左手系」に合わせて解釈しています。

*9:Unityと同じく、ZXYの順番で回転させています（参考）

*10:詳しい話は Klein(2012), "On the ellipsoid and plane intersection equation"

*11:実際試したわけではないです。

*12:実際は \sin^2 2\thetaの値