update WCS info

目的

\(N\) 天体のピクセル座標と天体座標のペアから画像の WCS 情報を更新することを考えます. ここでは以下のことを仮定します.

マッチングによって \(N\) 分の座標ペアを持っている,
チップあたりの歪み補正は必要ない,
ピクセルスケールは正しく求まっている.

天体座標とピクセル座標の対応

天体カタログ等に記載されている座標 \((\alpha,\delta)\) は fits 画像のヘッダ情報によって平面投影座標 \((\xi,\eta)\) へと変換されます. この変換には WCSlib などのライブラリを使用することができます.

\[ \left(\begin{array}{c} \alpha \\ \delta \end{array}\right) \xleftarrow{\text{WCSlib}} \left(\begin{array}{c} \phi \\ \theta \end{array}\right) \xleftarrow{\text{WCSlib}} \left(\begin{array}{c} \xi \\ \eta \end{array}\right) \xleftarrow{\text{affine}} \left(\begin{array}{c} x \\ y \end{array}\right) \]

天体の平面投影座標と画像上で検出されたソースの重心位置(ピクセル座標) \((x,y)\) はアフィン変換で結ばれます.

ピクセルスケールが既知のケース

まずはピクセルスケールが正しく求まっている場合を考えます. \((\xi,\eta)\) と \((x,y)\) の間には

\[ \left(\begin{array}{c} \xi \\ \eta \end{array}\right) = \left(\begin{array}{cc} \text{CD}_{1} & 0 \\ 0 & \text{CD}_{2} \end{array}\right) \left(\begin{array}{cc} \text{PC}_{11} & \text{PC}_{12} \\ \text{PC}_{21} & \text{PC}_{22} \end{array}\right) \left(\begin{array}{c} x - \text{CRPIX}_1\\ y - \text{CRPIX}_2 \end{array}\right) \]

という関係が成り立ちます. 簡単のために \(\text{CD}_n = 1\) となるように規格化します. \(n\) 番目の天体の座標ベクトルを \(\vec{\xi}_n\) および \(\vec{x}_n\) として

\[ \chi^2 \equiv \sum_n\left| \vec{\xi}_n - R_{\theta}\vec{x}_n - T \right|^2, \]

という量を最小化することが目的になります. ここで \(R_{\theta}\) および \(T\) はそれぞれ \(\theta\)-回転行列と平行移動成分を表しています. ここでは直接 \(\theta\) について最小化するのではなく, 回転行列を以下のように定義し直します.

\[ \begin{pmatrix} \xi \\ \eta \end{pmatrix} = \begin{pmatrix} c & {-}s \\ s & c \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} + \begin{pmatrix} T_x \\ T_y \end{pmatrix} \]

\(c^2 + s^2 = 1\) という制約の下で残差の二乗和 (\(\chi^2\)) を最小化することを考えます. ラグランジュの未定係数 \(\lambda\) を導入すると以下のように書けます.

\[ \begin{aligned} \chi^2 =& \sum_n\left|\xi_n-cx_n+sy_n-T_x\right|^2 \\ & + \sum_n\left|\eta_n-sx_n-cy_n-T_y\right|^2 + \lambda\left(c^2+s^2-1\right). \end{aligned} \]

これを \(c\), \(s\), \(T_x\), \(T_y\), \(\lambda\) について微分して最小値を求めます.

表記を簡素にするために以下の量を導入します.

\[ \begin{gathered} \overline{\xi x} = \sum_n (\xi_n x_n), \quad \overline{\eta x} = \sum_n (\eta_n x_n), \quad \overline{\xi y} = \sum_n (\xi_n y_n), \quad \overline{\eta y} = \sum_n (\eta_n y_n), \\ \overline{x x} = \sum_n (x_n x_n), \quad \overline{y y} = \sum_n (y_n y_n), \\ \overline{\xi} = \sum_n (\xi_n), \quad \overline{\eta} = \sum_n (\eta_n), \quad \overline{x} = \sum_n (x_n), \quad \overline{y} = \sum_n (y_n). \end{gathered} \]

微分係数はそれぞれ以下のように書けます.

\[ \def\xxbar{\overline{x x}} \def\xybar{\overline{x y}} \def\yybar{\overline{y y}} \def\xixbar{\overline{\xi x}} \def\etaxbar{\overline{\eta x}} \def\xiybar{\overline{\xi x}} \def\etaybar{\overline{\eta y}} \def\xbar{\overline{x}} \def\ybar{\overline{y}} \def\xibar{\overline{\xi}} \def\etabar{\overline{\eta}} \def\dd#1{\frac{\partial{\chi^2}}{\partial{#1}}} \begin{aligned} \frac{1}{2}\dd{c} &= (\xxbar+\yybar+\lambda)c + \xbar{T_x} + \ybar{T_y} - \xixbar - \etaybar, \\ \frac{1}{2}\dd{s} &= (\xxbar+\yybar+\lambda)s - \ybar{T_x} + \xbar{T_y} - \etaxbar + \xiybar, \\ \frac{1}{2}\dd{T_x} &= -\xibar + \xbar{c} - \ybar{s} + NT_x, \\ \frac{1}{2}\dd{T_y} &= -\etabar + \xbar{s} + \ybar{s} +NT_y, \\ \dd{\lambda} &= c^2 + s^2 -1. \end{aligned} \]

微分係数をすべて 0 とおいて \(T_x\), \(T_y\) を削除します.

\[ \def\xxbar{\overline{x x}} \def\xybar{\overline{x y}} \def\yybar{\overline{y y}} \def\xixbar{\overline{\xi x}} \def\etaxbar{\overline{\eta x}} \def\xiybar{\overline{\xi y}} \def\etaybar{\overline{\eta y}} \def\xbar{\overline{x}} \def\ybar{\overline{y}} \def\xibar{\overline{\xi}} \def\etabar{\overline{\eta}} \begin{aligned} \left[N(\xxbar+\yybar+\lambda)-(\xbar\,\xbar+\ybar\,\ybar)\right]c &= N(\xixbar+\etaybar) - \xibar\,\xbar - \etabar\,\ybar, \\ \left[N(\xxbar+\yybar+\lambda)-(\xbar\,\xbar+\ybar\,\ybar)\right]s &= N(-\xiybar+\etaxbar) - \etabar\,\xbar + \xibar\,\ybar. \end{aligned} \]

ここで \(\Delta^2\) を以下のように定義します.

\[ \Delta^2 \equiv \left[ N( \overline{\xi x} {+} \overline{\eta y}) - \overline{\xi}\,\overline{x} - \overline{\eta}\,\overline{y} \right]^2 + \left[ N(-\overline{\xi y} {+} \overline{\eta x}) - \overline{\eta}\,\overline{x} + \overline{\xi}\,\overline{y} \right]^2. \]

すると \(c^2 + s^2 - 1 = 0\) より以下の式を得ます.

\[ \begin{multline} \left[N(\xxbar+\yybar+\lambda)-(\xbar\,\xbar+\ybar\,\ybar)\right]^2 (c^2 + s^2 -1) = 0 \\ \qquad\qquad\qquad\iff \Delta^2 = \left[N(\xxbar+\yybar+\lambda)-(\xbar\,\xbar+\ybar\,\ybar)\right]^2 \end{multline} \]

これを \(\lambda\) について解くと最終的に以下の式を得ます.

\[ \lambda = \frac{\overline{x}\,\overline{x}{+}\overline{y}\,\overline{y} - N(\overline{x x} {+} \overline{y y}) \pm \Delta}{N}. \]

この \(\lambda\) を利用すると \(c\), \(s\), \(T_x\), \(T_y\) の解を以下の標式で得ることができます.

\[ \begin{aligned} &c = \frac{\overline{\xi}\,\overline{x} + \overline{\eta}\,\overline{y} - N( \overline{\xi x} {+} \overline{\eta y})}{\overline{x}\,\overline{x} + \overline{y}\,\overline{y} - N(\overline{x x} {+} \overline{y y} + \lambda)}, \\ &s = \frac{-\overline{\xi}\,\overline{y} + \overline{\eta}\,\overline{x} - N(-\overline{\xi y} {+} \overline{\eta x})}{\overline{x}\,\overline{x} + \overline{y}\,\overline{y} - N(\overline{x x} {+} \overline{y y} + \lambda)}, \\ &T_x = \frac{\overline{\xi} - c \overline{x} + s \overline{y}}{N}, \\ &T_y = \frac{\overline{\eta} + s \overline{x} + c \overline{y}}{N}. \end{aligned} \]

WCS のパラメタ PCi_j と CRPIXi は以下の式を満たします.

\[ \begin{pmatrix} \xi/\text{CDELT1} \\ \eta/\text{CDELT2} \end{pmatrix} = \begin{pmatrix} \text{PC1_1} & \text{PC1_2} \\ \text{PC2_1} & \text{PC2_2} \end{pmatrix} \begin{pmatrix} x - \text{CRPIX1}\\ y - \text{CRPIX2} \end{pmatrix}. \]

よって WCS パラメタを以下のように更新することができます.

\[ \begin{align} \text{PC1_1} &= c \\ \text{PC1_2} &= -s \\ \text{PC2_1} &= s \\ \text{PC2_2} &= c \\ \text{CRPIX1} &= - cT_x - sT_y \\ \text{CRPIX2} &= + sT_y - cT_y \\ \end{align} \]

ピクセルスケールが未知のケース

次にピクセルスケールも同時に最適化する場合を考えます. \((\xi,\eta)\) と \((x,y)\) には以下の関係が成り立ちます.

\[ \left(\begin{array}{c} \xi \\ \eta \end{array}\right) = \left(\begin{array}{cc} \text{CD}_{1} & 0 \\ 0 & \text{CD}_{2} \end{array}\right) \left(\begin{array}{cc} \text{PC}_{11} & \text{PC}_{12} \\ \text{PC}_{21} & \text{PC}_{22} \end{array}\right) \left(\begin{array}{c} x - \text{CRPIX}_1\\ y - \text{CRPIX}_2 \end{array}\right) \]

先ほどのケースと同様に残差の二乗和を最小化します. ここで \(R'_{\theta}\) はピクセルスケールの変換を含めた回転行列とします.

\[ \chi^2 \equiv \sum_n\left| \vec{\xi}_n - R'_{\theta}\vec{x}_n - T \right|^2, \]

今回は変換行列を以下ように定義します.

\[ \left(\begin{array}{c} \xi \\ \eta \end{array}\right) = \left(\begin{array}{cc} c_{11} & c_{12} \\ c_{21} & c_{22} \end{array}\right) \left(\begin{array}{c} x \\ y \end{array}\right) + \left(\begin{array}{c} T_x \\ T_y \end{array}\right) \]

\(c_{11}c_{21} + c_{12}c_{22} = 0\) という制約の下で \(\chi^2\) を最小化します. 再びラグランジュの未定係数 \(\lambda\) を導入します.

\[ \begin{aligned} \chi^2 =& \sum_n\left|\xi_n-c_{11}x_n-c_{12}y_n-T_x\right|^2 \\ & + \sum_n\left|\eta_n-c_{21}x_n-c_{22}y_n-T_y\right|^2 + \lambda\left(c_{11}c_{21}+c_{12}c_{22}\right). \end{aligned} \]

解析的に解こうとすると \(\lambda\) の 6 次式になったので, 解析的に解くことはあきらめて iterative に解くことにします. 次の関数を定義します.

\[ [x,y] = \overline{xy} - \frac{\overline{x}~\overline{y}}{N}. \]

微分係数が 0 になるという条件から \(\vec{c} = [c_{11}, c_{12}, c_{21}, c_{22}]^{t}\) は次の式を満たします.

\[ \begin{pmatrix} [x,x] & [x,y] & 0 & 0 \\ [y,x] & [y,y] & 0 & 0 \\ 0 & 0 & [x,x] & [x,y] \\ 0 & 0 & [y,x] & [y,y] \end{pmatrix} \vec{c} = \begin{pmatrix} [\xi,x] \\ [\xi,y] \\ [\eta,x] \\ [\eta,y] \end{pmatrix}. \]

この式を \(A\vec{c} = \vec{b}\) とおきます. 最適解は \(\vec{c}_0 = A^{-1}\vec{b}\) で得られます.

しかしこの最適解は拘束条件 \(c_{11}c_{21} + c_{12}c_{22} = 0\) を満たしているとは限りません. そこで \(\vec{c}\) の更新を以下の式で定義します.

\[ \vec{c}_{n+1} = (1-\varepsilon)\vec{c}_{n} + \varepsilon\vec{c}_\star. \]

\(\varepsilon > 0\) は更新の速度を決める定数です. 初期条件である \(\vec{c}_0\) は拘束条件を満たすように決めておきます.

\(f(\vec{c}) = c_{11}c_{21} + c_{12}c_{22}\) とおくと,

\[ \begin{aligned} {\rm d}f(\vec{c}) & = {\rm d}c_{11}{\cdot}c_{21} + {\rm d}c_{12}{\cdot}c_{22} + {\rm d}c_{21}{\cdot}c_{11} + {\rm d}c_{22}{\cdot}c_{12}, \\ & = [c_{21}, c_{22}, c_{11}, c_{12}] {\cdot} \vec{c}_\star, \\ & = \vec{c}^\dagger_n \cdot \vec{c}_\star. \end{aligned} \]

\(\vec{c}_\star\) を \(\vec{c}^\dagger_n\) と直交するように定めると \(f(\vec{c}_n) = 0 ~(\forall n \ge 0)\) を満たします. この更新方法で最適解を求めれば, 拘束条件の範囲内で解を探索することができます. \(\vec{c}_\star\) を以下の式で定義します.

\[ \vec{c}_\star = \vec{c}_0 - \frac{(\vec{c}^\dagger_n \cdot \vec{c}_0)} {|\vec{c}^\dagger_n|^2}\vec{c}^\dagger_n. \]

iterative に \(\vec{c}_n\) を更新して収束したところで解として採用します. \(\vec{c}\) が定まれば \(T_x, T_y\) は以下の式で得られます.

\[ T_x =\frac{\overline{\xi} - c_{11} \overline{x} - c_{12} \overline{y}}{N}, \qquad T_y =\frac{\overline{\eta} - c_{21} \overline{x} - c_{22} \overline{y}}{N}, \]

検算用コード

# fR: 回転行列用関数
fR = @(t) [cos(t) -sin(t); sin(t) cos(t)];

# D: スケール, R: 回転行列, T: 並進移動
D = diag(1+.1*randn(1,2));
R = fR(pi*rand());
T = 5*randn(1,2);

# Xi: インプット, Xt: 変換後の座標
Xi = 100*rand(200,2);
Xt = bsxfun(@plus, ((D*R)*Xi')'+.2*randn(200,2), T);

# 一次方程式 Ac = b を解く
a = @(Xi) [f(Xi(:,1),Xi(:,1)), f(Xi(:,1),Xi(:,2)); ...
           f(Xi(:,2),Xi(:,1)), f(Xi(:,2),Xi(:,2))];
A = @(X) blkdiag (a(X), a(X));
b = @(Xt,Xi) [f(Xt(:,1),Xi(:,1)); f(Xt(:,1),Xi(:,2)); ...
              f(Xt(:,2),Xi(:,1)); f(Xt(:,2),Xi(:,2))];

# 拘束条件を満たしつつ更新する
ct = @(c) [c(3) c(4) c(1) c(2)];
p = @(c,x) x-c*(ct(c)*(x))/(sumsq(c));
update = @(Xt,Xi,c,ep) (1.0-ep)*c + ep*p(c,A(Xi)\b(Xt,Xi));

# 初期条件を与えて iterative に解く
c = [1;0;0;1]; for n=2:100;
  c(:,n) = update(Xt,Xi,c(:,n-1),0.2);
  if(sumsq(diff(c(:,[n-1 n])')')<eps()) break; endif;
endfor

# Tx, Ty を解く関数
tx = @(Xt,Xi,c) [sum(Xt(:,1))-c(1)*sum(Xi(:,1)) ...
                -c(2)*sum(Xi(:,2))]/numel(Xt(:,1));
ty = @(Xt,Xi,c) [sum(Xt(:,2))-c(3)*sum(Xi(:,1)) ...
                -c(4)*sum(Xi(:,2))]/numel(Xt(:,2));

# 一次変換を適用する
mv = @(Xt,Xi,c) ...
   bsxfun(@plus, ([c(1) c(2); c(3) c(4)]*(Xi'))', ...
   [tx(Xt,Xi,c) ty(Xt,Xi,c)]);

# 結果を表示
plot(Xi(:,1),Xi(:,2),'o', ...
     Xt(:,1),Xt(:,2),'x', ...
   mv(Xt,Xi,c0)(:,1),mv(Xt,Xi,c0)(:,2),'+');

更新した WCS パラメタは以下の式で与えられます.

\[ \begin{aligned} \text{CDELT1} &= {-}\sqrt{c_{11}^2 + c_{12}^2}, \\ \text{CDELT2} &= {+}\sqrt{c_{21}^2 + c_{22}^2}, \\ \text{CRPIX1} &= -\left({+}c_{22} T_x - c_{12} T_y\right) \mathbin{/} \det{C}, \\ \text{CRPIX2} &= -\left({-}c_{21} T_y + c_{11} T_y\right) \mathbin{/} \det{C}, \\ \text{PC1_1} &= c_{11}\mathbin{/}\text{CDELT1}, \\ \text{PC1_2} &= c_{12}\mathbin{/}\text{CDELT1}, \\ \text{PC2_1} &= c_{21}\mathbin{/}\text{CDELT2}, \\ \text{PC2_2} &= c_{22}\mathbin{/}\text{CDELT2}, \end{aligned} \]