実装して学ぶデータ処理システム その4 - 座標変換の実際

構造 プログラミング

前回(本記事公開時には未公開)、Kabsch 変換の意味について説明した。今回は、実装上の問題点について述べる。DIALS のレポジトリにある James の文章 が参考になる。

逆空間の一点とKabsch 変換後の反射座標系での一点は、数学的には一対一に対応している。しかし、プログラムとして実装する上では、それぞれを有限の大きさを持ったボクセルの集合として表現するので、多対多の関係になってしまう。つまり、反射座標系の 1 つのボクセルに対して、複数のフレーム上の複数のピクセルが対応しうる。逆に、あるフレームのあるピクセルが、反射座標系では複数のボクセルにまたがることもある。

このような複雑な対応関係を処理するために、XDS は複数の戦略を取ってきた。Kabsch, Wolfgang. "Evaluation of single-crystal X-ray diffraction data from a position-sensitive detector" Journal of Applied Crystallography 21.6 (1988) で述べられている初期のバージョンでは、1 つの検出器ピクセルの値は、反射座標系で近傍の 2 * 2 * 2 = 8 ボクセルに線形補間で配られる。たとえば、検出器ピクセルの中心座標が、プロファイル座標系で (3.2, 2.4, 1.5) に対応するならば、(3, 2, 1) から (4, 3, 2) までの 8 つのボクセルに配るのである。この方法は、反射座標系の grid の刻みが細かくなって、1つの検出器ピクセルがもっと多くのボクセルにまたがるようになると破綻してしまう。

2010 年の論文 Kabsch, Wolfgang. "Integration, scaling, space-group assignment and post-refinement" Acta Crystallographica Section D (2010): 133-144 で説明されているように、新しいバージョンの XDS では、1 つの検出器ピクセルを 5 * 5 = 25 個のサブピクセルに分割して、それぞれを対応するボクセルへ配るようになった。ただし、各ボクセルの間で補間することはしない。この方法でも、1 つの検出器ピクセルが 25 個以上のボクセルにまたがるとうまくいかないはずだが、実用上は問題ないようである(BEAM_DIVERGENCE などをうまく取っている?)。

理論的には、1 つの検出器ピクセルは、反射座標系で凸多面体に対応する。この多面体と交差するボクセルを列挙し、交差する体積を解析的に求めることは、複雑だが可能である。DIALS における XDS アルゴリズムの実装は、この方向で進められていると聞いている。＃対応するコードと資料を要確認

XDS では、プロファイル座標系でのボクセル数は NUMBER_OF_PROFILE_GRID_POINTS_ALONG_ALPHA/BETA, NUMBER_OF_PROFILE_GRID_POINTS_ALONG_GAMMA パラメータで設定される。デフォルトは 9 で、最大は 21 だ。これがプロファイル空間でどれだけの幅にマップされるかを決めるのが BEAM_DIVERGENCE と REFLECTING_RANGE である。 BEAM_DIVERGENCE_E.S.D と REFLECTING_RANGE_E.S.D は、三次元正規分布で表される解析的なプロファイルの幅(標準偏差)であり、強い反射のスポット形状から決定される。両側に 3 sigma 取れば正規分布の 99 % をカバーするから、CUT が 2 の場合、 BEAM_DIVERGENCE と REFLECTING_RANGE は 各 E.S.D. の 6 倍あれば十分である。実際には、正規分布よりも広がりがある(over dispersion) かもしれないから、少し多めにとっておいたほうが無難である。実際の XDS がどのような基準で決めているかは論文では明記されていない(邪悪!)が、実は 7 倍らしい。

実装して学ぶデータ処理システム その3 - Kabsch transform の意味

構造 プログラミング

今回は XDS における最大の特徴、Kabsch transform の意味について説明する。論文の 2.3 節に対応する。

Kabsch transform は、goniostat 座標系から reflection 座標系への変換である。これまでに出てきた検出器系や実験室系との変換と違い、この変換は局所的である。局所的というのは、1つ1つの逆格子点ごとに固有の reflection 座標系があって、その点の近傍でしか変換が意味をなさないということだ。数学が得意な人は、多様体における局所座標系の貼りあわせをイメージするかもしれない。ただ、検出器上の1つのピクセルは、最近傍の反射がもつ局所座標系だけに属することになっており(pixel labeling)、局所座標系同士が重なりあうことはないから、多様体よりも単純である。

ある反射 hkl の固有座標系は、逆空間(goniostat 座標系で考えるのがよい)でその反射に対応する格子点を原点とし、基底ベクトル

\mathbf{e_1} = \frac{\mathbf{S} \times \mathbf{S_0}}{|\mathbf{S} \times \mathbf{S_0}|}
\mathbf{e_2} = \frac{\mathbf{S} \times \mathbf{e_1}}{|\mathbf{S} \times \mathbf{e_1}|}
\mathbf{e_3} = \frac{\mathbf{S} + \mathbf{S_0}}{|\mathbf{S} + \mathbf{S_0}|}

によって定義される。1つ目の式は、\mathbf{e_1} は \mathbf{S}, \mathbf{S_0} に直交する長さ1のベクトルだと言っている。二つ目の式は 、こうしてできた \mathbf{e_1} と \mathbf{S} に直交する長さ1のベクトルが \mathbf{e_2} だと言っている。どちらも散乱ベクトル \mathbf{S} に直交するから、この2つのベクトルは Ewald 球の接平面上にある。\mathbf{e_3} の向きを考えるために、\mathbf{e_1} との内積を取ると、

\mathbf{e_1} \cdot \mathbf{e_3} = \frac{(\mathbf{S} \times \mathbf{S_0}) \cdot (\mathbf{S} + \mathbf{S_0})}{|\mathbf{S} \times \mathbf{S_0}||\mathbf{S} + \mathbf{S_0}|} = \frac{(\mathbf{S} \times \mathbf{S_0}) \cdot \mathbf{S} + (\mathbf{S} \times \mathbf{S_0}) \cdot \mathbf{S_0}}{|\mathbf{S} \times \mathbf{S_0}||\mathbf{S} + \mathbf{S_0}|} = 0

となるから、\mathbf{e_1} と \mathbf{e_3} は直交している。同様に \mathbf{p^*} = \mathbf{S} - \mathbf{S_0} との内積をとっても、

\mathbf{p^*} \cdot \mathbf{e_3} = \frac{(\mathbf{S} - \mathbf{S_0}) \cdot (\mathbf{S} + \mathbf{S_0})}{|\mathbf{S} + \mathbf{S_0}|} = \frac{\mathbf{S} \cdot \mathbf{S} + \mathbf{S} \cdot \mathbf{S_0} - \mathbf{S_0} \cdot \mathbf{S} - \mathbf{S_0} \cdot \mathbf{S_0}}{|\mathbf{S} + \mathbf{S_0}|} = 0

となって、こちらも直交している。ここで、内積が可換であることと、|\mathbf{S}| = |\mathbf{S_0}| (Laue 条件) を用いた。

反射座標系の成分のスケールとしては、角度(度)を用いることになっている。長さの単位として角度というのは不思議に感じられるかもしれない。これが可能となるのは、\mathbf{e_1} と \mathbf{e_2} は Ewald 球の接平面上にあり、ごく近傍では Ewald 球の表面(曲面)とも同一視できるからである。Ewald 球面上の弧の長さは、半径 |\mathbf{S_0}| = \frac{1}{\lambda}と中心角(ラジアン)をかけたものだから、長さの基準として角度を用いても問題ないのだ。

数式はこのくらいにしておいて、意味について考えよう。検出器上で観察される spot profile は、still shot の場合、次のような影響をうける。

  • Ewald 球の曲率に由来する広がり (高角のほうが大きい)
  • スペクトル幅による広がり (高角のほうが大きい)
  • sensor 厚みによる位置のズレ

振動写真の場合は、さらに次のような影響も加わる。

  • 回転軸から遠い( \rho が大きい)反射は、大きい線速度を持つ

これらを理解するために、まず、逆空間における1つ1つの「反射球」は、一定の大きさを持っていることを思い出そう。1つの単位胞をフーリエ変換したもの、つまり 連続的な molecular transform に無限サイズの格子に由来する逆格子を掛け算したものが理想結晶の逆空間を与える。実際の結晶は有限のサイズを持つから、その shape transform を逆格子点の各点に畳み込んだものとなる。詳しくは◯◯を参照。

一方、検出器上で観測される現象は、#未稿