CrystFEL のピークサーチについて

CrystFEL におけるピークサーチについてソースコードを読んだメモ。このエントリはGPLとする。

実際のワークフローでは、Cheetah などのより高速なソフトウェアで画像に反射が含まれる(hit) かどうかを判定し、hit だけを CrystFEL に食わせることになる。

libcrystfel/src/peaks.c にコードがある。入り口は

void search_peaks(struct image *image, float threshold, float min_gradient,
                  float min_snr, double ir_inn, double ir_mid, double ir_out,
                  int use_saturated)

だ。image は複数のパネルからなるが、パネルをまたがるような反射は無視することになっているので、パネルごとに search_peaks_in_panel を呼び出して処理している。処理後のスポットリストは、image->features に記入される。

static void search_peaks_in_panel(struct image *image, float threshold,
                                  float min_gradient, float min_snr,
                                  struct panel *p,
                                  double ir_inn, double ir_mid, double ir_out,
                                  int use_saturated)

は、パネル内のスポットを処理する。1つのパネルに対応する領域は、アセンブリ後のイメージ image->data の一部分(（min_ss, min_fs)-(max_ss, max_fs) の範囲の矩形)をなす。その範囲の各ピクセルについて、四方の差から gradient を求め、閾値 min_gradient と比較する。閾値以上ならば、そのピクセルの周辺での極大を求める。極大が、開始点よりも ir_inn 以上離れていれば reject する。「周辺」と書いたが、その条件が少しおかしなループになっていて、効率も悪そうに見える。意図的なものかは不明。

次に、integrate_peak 関数を呼び出して、ピークの積分と重心の計算を行う。I/sigma が min_snr 未満なら reject。2 * ir_inn の範囲に他のスポットがあっても reject する。

int integrate_peak(struct image *image, int cfs, int css,
                   double *pfs, double *pss,
                   double *intensity, double *sigma,
                   double ir_inn, double ir_mid, double ir_out,
                   int *bgPkMask, int *saturated)

は、アセンブリ後のイメージ image における座標 (css, cfs) を含むピークの積分を行い、強度 intensity, 標準偏差 sigma, 重心 (pss, pfs) を返す。(p_css, p_cfs) は、パネル座標系。

まず、背景を見積る。(css, cfs) を中心として、距離が ir_mid と ir_out の間であるような環状領域をバックグラウンドと見なして、平均 bg_mean と分散 bg_var を計算する。次に、距離が ir_min 以下であるような円をピークと考えて、ピクセル値から上のバックグラウンド平均を引いたものを積分する。重心も求める。標準偏差 sigma は、背景の分散 bg_var と波長に依存する検出器ノイズ(?) aduph をピクセル分加算したものの平方根を取って求めている。

なお、cull_peaks(と、本体である cull_peaks_in_panel) は、画像の row や column に一列に連なるスポットを除去する後処理である。CCD の電荷漏れによる smear を除去するためのものかと思われるが、「使うな」と書いてあるので省略。