ジョーンズ計算法によるセナルモン法の説明

複屈折媒質(リターダ)に光が入射すると,直交する偏光成分間に位相差が生じる.この位相差のことをリタデーションといい,実験的に測定する手法としてセナルモン法 (Senarmont method) がある.意外とセナルモン法について書かれた書籍が少ないので,備忘録を兼ねてジョーンズ計算法 (Jones calculus) によりセナルモン法を説明してみる.

セナルモン法の手順

セナルモン法の手順は以下の通りである.割と面倒.

  1. レーザーの光軸上に偏光子を置く
  2. 光が透過しないように偏光子の回転角を調整する(偏光面と透過軸を直交させる)
  3. 偏光子の手前に {\lambda/4} 板を置き,光が透過しないように回転角を調整して(高速軸を偏光面と一致させ)一旦取り出す
  4. リターダを偏光板の手前に置き,光が透過しないように回転角を調整する(高速軸を偏光面と一致させる)
  5. リターダを {45^{\circ}} だけ回転させる
  6. リターダと偏光子の間に一旦取り出した {\lambda/4} 板を再び入れる
  7. 光が透過しないように偏光子の回転角を調整する.その回転角の2倍がリタデーションである

ジョーンズ行列による偏光素子の表記

偏光が絡む問題を数学的に扱う方法としてはジョーンズ計算法ミュラー計算法 (Mueller calculus) があるが,完全偏光だけを扱うのであればジョーンズ計算法が簡単でよい.

リターダのジョーンズ行列

リターダの高速軸と低速軸をそれぞれ {x, y} 軸として,リタデーション {\phi} を生じるリターダのジョーンズ行列 {\mathbf{J}(\phi)} は次の式で表せる.


  \begin{equation}
    \mathbf{J}(\phi) =
    \begin{bmatrix}
      \mathrm{e}^{\mathrm{i}\frac{\phi}{2}} & 0 \\
      0 & \mathrm{e}^{-\mathrm{i}\frac{\phi}{2}}
    \end{bmatrix}
  \end{equation}

高速軸が {x} 軸から {\theta} 傾いたリターダのジョーンズ行列 {\mathbf{J}(\phi, \theta)} は,2次元の回転行列{\mathbf{R}(\theta)} として {\mathbf{J}(\phi, \theta) = \mathbf{R}(-\theta)\mathbf{J}(\phi)\mathbf{R}(\theta)} で表せるので,


  \begin{equation}
    \mathbf{J}(\phi, \theta) = 
    \begin{bmatrix}
      \cos\dfrac{\phi}{2} + \mathrm{i}\sin\dfrac{\phi}{2}\cos 2\theta & \mathrm{i}\sin\dfrac{\phi}{2}\sin 2\theta \\
      \mathrm{i}\sin\dfrac{\phi}{2}\sin 2\theta & \cos\dfrac{\phi}{2} - \mathrm{i}\sin\dfrac{\phi}{2}\cos 2\theta
    \end{bmatrix}
  \end{equation}

となる.

{\lambda/4}長板と偏光子のジョーンズ行列

セナルモン法で用いる {\lambda/4} 板と偏光子のジョーンズ行列を以下に示しておく.


  \begin{align}
    \mathbf{J}_{\mathrm{QWP}}(\theta) = \mathbf{J}\left(\dfrac{\pi}{2}, \theta\right) &= 
    \begin{bmatrix}
      \dfrac{\sqrt{2}}{2} + \mathrm{i}\dfrac{\sqrt{2}}{2}\cos 2\theta & \mathrm{i}\dfrac{\sqrt{2}}{2}\sin 2\theta \\
      \mathrm{i}\dfrac{\sqrt{2}}{2}\sin 2\theta & \dfrac{\sqrt{2}}{2} - \mathrm{i}\dfrac{\sqrt{2}}{2}\cos 2\theta
    \end{bmatrix}
    \\
    \mathbf{J}_{\mathrm{POL}}(\theta) &=
    \begin{bmatrix}
      \cos^{2}\theta & \sin\theta\cos\theta \\
      \sin\theta\cos\theta & \sin^{2}\theta
    \end{bmatrix}
  \end{align}

偏光子の {\theta} は透過軸の {x} 軸からの傾きである.

ジョーンズ計算法によるセナルモン法の説明

入射光は {[1~0]^{T}} 偏光であるとする.上に書いたセナルモン法の手順は結局のところ,{[1~0]^{T}} 偏光の光をリターダ {\mathbf{J}(\phi, 45^{\circ})}{\lambda/4}{\mathbf{J}_{\mathrm{QWP}}(0^{\circ})},偏光子 {\mathbf{J}_{\mathrm{POL}}(\theta+90^{\circ})} の順に入射させることに相当するので,透過光 {\mathbf{u}} はジョーンズ計算法により,


  \begin{align}
    \mathbf{u}&=\mathbf{J}_{\mathrm{POL}}(\theta+90^{\circ})\mathbf{J}_{\mathrm{QWP}}(0^{\circ})\mathbf{J}(\phi, 45^{\circ})\begin{bmatrix}1 \\ 0\end{bmatrix}\\
    &=
    \begin{bmatrix}
      \sin^{2}\theta & -\sin\theta\cos\theta \\
      -\sin\theta\cos\theta & \cos^{2}\theta
    \end{bmatrix}
    \begin{bmatrix}
      \mathrm{e}^{\mathrm{i}\frac{\pi}{4}} & 0\\
      0 & \mathrm{e}^{-\mathrm{i}\frac{\pi}{4}}
    \end{bmatrix}
    \begin{bmatrix}
      \cos\dfrac{\phi}{2} & \mathrm{i}\sin\dfrac{\phi}{2} \\
      \mathrm{i}\sin\dfrac{\phi}{2} & \cos\dfrac{\phi}{2}
    \end{bmatrix}
    \begin{bmatrix}1 \\ 0\end{bmatrix}\\
    &=
    \begin{bmatrix}
      \mathrm{e}^{\mathrm{i}\frac{\pi}{4}}\sin\theta\left(\cos\dfrac{\phi}{2}\sin\theta-\sin\dfrac{\phi}{2}\cos\theta\right) \\
      \mathrm{e}^{\mathrm{i}\frac{\pi}{4}}\cos\theta\left(\sin\dfrac{\phi}{2}\cos\theta-\cos\dfrac{\phi}{2}\sin\theta\right)
    \end{bmatrix}\\
    &=
    \begin{bmatrix}
      \mathrm{e}^{\mathrm{i}\frac{\pi}{4}}\sin\theta\sin\left(\theta-\dfrac{\phi}{2}\right) \\
      -\mathrm{e}^{\mathrm{i}\frac{\pi}{4}}\cos\theta\sin\left(\theta-\dfrac{\phi}{2}\right)
    \end{bmatrix}
  \end{align}

と表され,光強度は


  \begin{align}
    |\mathbf{u}|^{2}=\sin^{2}\left(\theta-\dfrac{\phi}{2}\right)
  \end{align}

となる.上の手順 7. は {\sin^{2}(\theta-\phi/2)=0} となる {\theta} を見つけることに他ならない.これは当然 {\theta=\phi/2} なので,偏光子の回転角 {\theta} の2倍がリタデーション {\phi} となる.

シングルピクセルイメージング

イメージングというと普通は CCD や CMOS などのアレイセンサを使うが,フォトダイオードのように空間分化能を持たない点型検出器を使ったイメージング手法を圧縮イメージングとかシングルピクセルイメージング (SPI) という.この方法の何がうれしいかというと,高感度なので SNR が低くても解像可能らしい.他にも点型検出器を使うので普通のイメージセンサが感度を持たないテラヘルツ波などでイメージングできる,結像光学系が必要ない,といった利点もある.シングルピクセルイメージングの基本的なアイデアは空間的に変調された照明光を使う点にある.

問題の定式化

照明光の空間パターンを1次元に並べたベクトルを {\mathbf{a}_{i}},イメージング対象物の透過率/反射率の空間分布を1次元に並べたベクトルを {\mathbf{x}},透過率/反射率の光強度を {b_{i}} とすると,


  \begin{equation}
    \langle\mathbf{a}_{i}, \mathbf{x}\rangle = \mathbf{a}_{i}^{T}\mathbf{x} = b_{i}
  \end{equation}

となるのはいいだろう.あるいは連続系で表した方が分かりやすいかもしれない.


  \begin{equation}
    \iint a_{i}(\mathbf{r})x(\mathbf{r})\mathrm{d}^{2}r =b_{i}
  \end{equation}

添字の {i=1, 2, \dots} は測定回を表す.空間パターンの生成には空間光変調器 (SLM) やデジタルミラーデバイス (DMD) などが使われるようである.{\mathbf{a}_{i}^{T}} を列方向に並べた行列を {\mathbf{A}}{b_{i}} を並べたベクトルを \mathbf{b} とすると,まとめて


  \begin{equation}
    \mathbf{A}\mathbf{x} = \mathbf{b}
  \end{equation}

と書ける.何のことはない,線形方程式系である.ただしこの線形方程式系は普通は劣決定系である.というのも,対象物の解像度を {N^{2}},測定回数を {M} とすると,{\mathbf{A}}{M}N^{2} 列となるが,以下に説明するように {N^{2}>M} でもイメージング可能なのである.それに(優)決定系にしようと思えば,解像度を少し高くしただけで多数の測定が必要になることはすぐ分かる.

解析的アプローチ

劣決定系の線形方程式系の解法としてすぐ思い浮かぶのは,正規方程式から {\|\mathbf{A} \mathbf{x}-\mathbf{b}\|_{2}} を最小化する {\mathbf{x}} 求めるという方法である.この場合,普通 {\mathbf{A}^{T}\mathbf{A}}正則行列ではないので,最急降下法共役勾配法等の反復法が使われる.反復法で最小二乗解を求める方法はそこそこ上手くいくが,測定結果にノイズが乗っている場合いわゆる「過剰適合」を生じる.そこでもう少し凝った方法として正則化を使う手がある.スパースモデリングに基づく方法だと,{\mathbf{x}} がある種の基底変換でスパース(ほとんどがゼロ成分)になるという性質を利用する.基底変換の行列を {\mathbf{D}},変換後のベクトル(基底に関する成分)を {\mathbf{c}} とすると,{\mathbf{D}\mathbf{x}=\mathbf{c}} となるが,{\mathbf{c}}{\ell_{1}}-norm を最小化してやるのである(本当は {\ell_{0}}-norm を最小化すべきだが,この場合は組合せ最適化の問題となってしまう.なお,{\ell_{0}}-norm は通常の意味でのノルムではない).問題としては次の様に定式化できる.


  \begin{equation}
    \mathrm{min}~\|\mathbf{c}\|_{1}\quad\mathrm{s. t.}\quad\mathbf{A}\mathbf{x}=\mathbf{b},~\mathbf{D}\mathbf{x}=\mathbf{c}
  \end{equation}

この手の問題は Alternating Direction Method of Multipliers (ADMM) で割と効率的に解ける(「割と」と書いたのは逆行列を求める必要があるため).ADMM では次の拡張ラグランジュ関数の最小化問題を解く.


  \begin{equation}
    \mathrm{min}~L=\|\mathbf{c}\|_{1}+\mathbf{h}^{T}_{1}(\mathbf{D}\mathbf{x}-\mathbf{c})+\frac{\mu_{1}}{2}\|\mathbf{D}\mathbf{x}-\mathbf{c}\|_{2}+\mathbf{h}^{T}_{2}(\mathbf{A}\mathbf{x}-\mathbf{b})+\frac{\mu_{2}}{2}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_{2}
  \end{equation}

ここで {\mathbf{h}}ラグランジュ乗数,{\mu} は調整パラメータである.{\mathbf{D}} には {\mathbf{x}} の勾配を取る行列を使うとよい結果が得られるようである(この場合,{\|\mathbf{c}\|_{1}} が Total Variation となる).

統計光学的アプローチ

上記とは違ったアプローチとして {\mathbf{A}}\mathbf{b} の相関を利用する方法がある.詳しい原理はこちらを見てもらうとして,手順としては次の計算を行うだけである.


  \begin{equation}
    \mathbf{g} = \frac{1}{M}\mathbf{A}^{T}(\mathbf{b}-\langle\mathbf{b}\rangle)
  \end{equation}

ここで {\langle\rangle} はアンサンブル平均である.この計算は {\mathbf{a}_{i}} のゆらぎ {\Delta\mathbf{a}_{i}=\mathbf{a}_{i}-\langle\mathbf{a}_{i}\rangle}と,{b_{i}} のゆらぎ {\Delta b_{i}=b_{i}-\langle b_{i}\rangle}2次コヒーレンス {\langle\Delta\mathbf{a}_{i}\Delta b_{i}\rangle} を求めることに相当する.この方法は特にゴーストイメージングと呼ばれているようである.反復計算を必要としないのはよいが,{M\to\infty} での誤差の漸近挙動があまりよろしくない.ただしノイズには強い.

直交基底を使う方法

これまで係数行列(照明光の空間パターン){\mathbf{A}} について特に言及しなかったが,こちらを工夫する手もある.何も考えずランダムなパターンを用いるという手もあるが,直交基底を用いると空間周波数ごとに測定できて効率的である.実際,アダマール行列に基づいた空間パターンを使うことでゴーストイメージングより鮮明度の高い像が得られるようである.あるいは事後的に {\mathbf{A}} を直交化するという方法も考えられる.グラム・シュミットの方法で {\mathbf{A}} を正規直交化すると({\mathbf{b}} にも同じ操作を施す),かなりよい像が得られるようである.ただしこの方法ではノイズに弱くなる模様.

回折の影響

変調器で生成された照明光は回折するので,必要に応じてフレネル回折やら角スペクトル法で回折パターンを求める必要がある.ただしリレー光学系を組み込めばこの処理は必要ない.

参考文系

あとで書く.

ベクトル化した2次元データに対する離散フーリエ変換行列

機械学習では画像のような2次元のデータをベクトル化して1次元のデータとして扱うことが多い.このようにベクトル化した2次元データに対する離散フーリエ変換 (DFT) 行列を定義できれば便利である.

長さ {N} のベクトル {\mathbf{x}} に対する1次元の DFT は,{\mathbf{W}_{N}}{N{\times}N}DFT 行列として次の様に表せる.


  \begin{equation}
    \mathrm{DFT}_{1}(\mathbf{x}) = \mathbf{W}_{N}\mathbf{x}
  \end{equation}

同様に {M}{N} 列の2次元データ {\mathbf{X}} に対する2次元の DFT は次の様に表せる.


  \begin{equation}
    \mathrm{DFT}_{2}(\mathbf{X}) = \mathbf{W}_{M}\mathbf{X}\mathbf{W}_{N}
  \end{equation}

さて,2次元のデータ(行列)をベクトル化する操作を {\mathrm{vec}(\cdot)} と表すと,次の関係が成り立つ.


  \begin{equation}
    \mathrm{vec}(\mathbf{AXB}) = \left(\mathbf{B}^{T}\otimes\mathbf{A}\right)\mathrm{vec}(\mathbf{X})
  \end{equation}

ここで {\otimes}クロネッカー積である.2番目の式とあわせて考えると,次の関係を得る*1


  \begin{align}
    \mathrm{vec}(\mathrm{DFT}_{2}(\mathbf{X})) &= \mathrm{vec}(\mathbf{W}_{M}\mathbf{X}\mathbf{W}_{N})\\
    &= \left(\mathbf{W}_{N}^{T}\otimes\mathbf{W}_{M}\right)\mathrm{vec}(\mathbf{X})\\
    &= \left(\mathbf{W}_{N}\otimes\mathbf{W}_{M}\right)\mathrm{vec}(\mathbf{X})\\
  \end{align}

というわけで,ベクトル化した2次元データに対する DFT 行列は {\mathbf{W}_{N}\otimes\mathbf{W}_{M}} で定義できることが分かる.ちなみに {\mathbf{X}}{128{\times}128},DFT 行列の要素を std::complex< double > とすると,この DFT 行列は 4 GB のメモリを必要とする.

*1:DFT 行列の性質 {\mathbf{W}^T=\mathbf{W}} に注意.

シルベスターの方法で生成したアダマール行列からウォルシュ型アダマール行列を得る方法 (C++)

信号処理や統計学で使われるアダマール行列は,以下のようにシルベスターの方法で再帰的に生成するのが(おそらく)一般的.


  \begin{align}
    \mathbf{H}_{2} &=\left[
    \begin{array}{rr}
      1 & 1 \\\
      1 & -1
    \end{array}\right]\\\
    \mathbf{H}_{2^{k}} &= \mathbf{H}_{2}\otimes\mathbf{H}_{2^{k-1}}
  \end{align}

ここで \otimesクロネッカー積を表す.例として \mathbf{H}_{8} を示す.


  \begin{equation}
  \mathbf{H}_{8} =\left[
  \begin{array}{rrrrrrrr}
    1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\\
    1 & -1 & 1 & -1 & 1 & -1 & 1 & -1\\\
    1 & 1 & -1 & -1 & 1 & 1 & -1 & -1\\\
    1 & -1 & -1 & 1 & 1 & -1 & -1 & 1\\\
    1 & 1 & 1 & 1 & -1 & -1 & -1 & -1\\\
    1 & -1 & 1 & -1 & -1 & 1 & -1 & 1\\\
    1 & 1 & -1 & -1 & -1 & -1 & 1 & 1\\\
    1 & -1 & -1 & 1 & -1 & 1 & 1 & -1
  \end{array}\right]
  \end{equation}

さて,この方法で生成したアダマール行列について,各行の符号反転の回数(周波数)を小さい順に並べ替えた方が都合がいい場合がある.そのような行列をウォルシュ型アダマール行列と呼んだりするが,次のような方法で並べ替えることができる.

/**
 * @brief ビット逆転(32 bit 用)
 * @param[in] a 32 bit 符号なし整数
 * @return ビット逆転した 32 bit 符号なし整数
 */
unsigned int bitreverse( const unsigned int a ) {
    unsigned int x = a;
    x = ( ( ( x & 0xaaaaaaaa ) >> 1 ) | ( ( x & 0x55555555 ) << 1 ) );
    x = ( ( ( x & 0xcccccccc ) >> 2 ) | ( ( x & 0x33333333 ) << 2 ) );
    x = ( ( ( x & 0xf0f0f0f0 ) >> 4 ) | ( ( x & 0x0f0f0f0f ) << 4 ) );
    x = ( ( ( x & 0xff00ff00 ) >> 8 ) | ( ( x & 0x00ff00ff ) << 8 ) );
    return( ( x >> 16 ) | ( x << 16 ) );
}

/**
 * @brief 並び替えのための配列生成
 * @param[in] k 行列の次数 (2^k)
 * @param[out] v 並び替えのための配列
 */
void permutation( const unsigned int k,
                  std::vector< unsigned int > &v ) {
    const unsigned int N = 1 << k;
    v.resize(N);
    for ( unsigned int i = 0; i < N; i++ ) {
        unsigned int j = bitreverse(i) >> ( 31 - k );
        j = j ^ ( j >> 1 );
        j = ( j << ( 32 - k ) ) >> ( 32 - k );
        v[i] = j;
    }
}

たとえばpermutation( 3, v )とした場合はv = [0 4 6 2 3 7 5 1]となるので,この順に行を並べ替えれば8次のウォルシュ型アダマール行列を得られる.上のプログラムは当然ながらsizeof( unsigned int ) == 4という前提.

Julia 1.0 を Jupyter Notebook から使う

以下は2018年8月11日時点の情報です。

Anaconda3 がインストールされている前提。Julia を起動して]キーを入力して Pkg REPL-mode に入る。

julia> ]

Conda→IJulia の順でパッケージをインストール。Conda には#master master ブランチの最新版を取得するようにする(#masterを付けないとビルドでこける)。

(v1.0) pkg> add Conda#master
(v1.0) pkg> add IJulia
julia> using IJulia

Jupyter Notebook を起動して "New" をクリック、"Julia 1.0.0" が表示されていれば OK。

f:id:ystt:20180811140952j:plain

2018年9月5日追記

#masterを付けなくてもビルドが通るようになった模様.rmでインストール済みのパッケージを削除してから,addで入れ直し.

Anaconda のアップデート

いつも忘れるのでメモ.

> conda update conda
> conda update anaconda
> conda update --all

conda→Anaconda→パッケージの順番でアップデート.Windows 版を使っていれば conda.exe は C:\ProgramData\Anaconda3\Scripts にある.

TeX Live で LuaJITLaTeX を使う (Windows)

毎回忘れるのでメモ.

  1. kpsewhich fmtutil.cnf で fmtutil.cnf の path を探す。
  2. fmtutil.cnf をエディタで開き、luajitlatex を検索。該当行の行頭は #! でコメントアウトされているので、これを消して行を有効にする。
  3. fmtutil-sys --byfmt luajitlatex を実行。
  4. luajittex --fmt=luajitlatex.fmt foo.texlatex として使用可。

Windows では上記の設定を行うと luajitlatex.exe が使えるようになる.

> luajitlatex.exe foobar.tex
> luajitlatex.exe -cmdx ほげほげ.tex

日本語のファイル名を扱う場合は,-cmdxオプションを第1引数に指定する必要がある.TeXstudio から呼び出す場合は,

luajitlatex.exe -cmdx -synctex=1 -interaction=nonstopmode %.tex

とでもしておけばよい.