実装して学ぶデータ処理システム その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 倍らしい。