OPM (天体カタログ同定ソフト)

このページでは、位置座標と等級が求められている2つの天体カタログに対して、 座標の変換式を求めて同定も行う自作のプログラムを公開し、 使用法などを説明しています。たとえば、測光ソフトでどのピクセルに何等の星が あるかというリストをつくった後に、カタログの赤道座標と比較して 座標の変換をしたり、カタログの等級と比較したりできます。 また、３つのフィルターでとった画像の座標系を決めて、合成カラー画像を 作ったりできます。

座標の変換式を求める根本部分は、Tabur (2007, PASA, 24, 189)によって 提案されたOptimistic Pattern Matching (OPM)というアルゴリズムに 基づいています。OPMアルゴリズムは、IRAFのxyxymatchタスクなどが 採用している従来の方法よりも速いと期待されています。 コードはC言語で書かれており、適当なC言語のライブラリがあれば、 結果を図示したり、求めた座標の変換式をそのままFITS画像のWCSに 入力したりもできます。

本ソフトは、学術的な目的であれば、自由に利用できます (著作権情報はこちら)。 論文などで引用、あるいは謝辞に入れてくださる場合は、 こちらを参考にしてください。

目次

参考文献

V. Tabur, "Fast algorithms for matching CCD images to a stellar catalogue",
Publications of the Astronomical Society of Australia, Volume 24, Issue 4, pp. 189-198. (2007)
(プレプリントがarXiv:0710.3618v1[astro-ph]で入手できます)
(Tabur自身によるソフトのページがこちらにあります)

更新情報

松永典之

ダウンロード

以下のファイルをダウンロードして、適当なディレクトリに置いてください。 ファイルサイズは約1.4MBです。
OPM1_4.tgz

インストールの方法

Unix系の計算機環境でのインストール方法は以下の通りです。

  1. tar zxvf OPM1_4.tgz
  2. cd OPM
  3. make

以下にverboseモード、plotモード、cfitsioモードの３つの種類の インストール方法を説明しています。verboseモードは、opmプログラムを 実行したときにいろいろな情報を画面に表示するかどうかで、 デフォルトではverboseモードになっています。Makefileを編集してから makeを行うと非verboseモードにすることもできます。 一方、plotモードとcfitsioモードは、opmの実行には 直接関係のないオプションのプログラムをコンパイルします。 PGPLOTやcfitsioのライブラリがある場合は、これらのオプションを インストールすることで、いくつか補助的なプログラムを実行できるようになります。

verboseモードのインストール

opm.cは、いろいろな情報をモニターに表示させたり、させなかったりすることが できます。Makefileの最初の方で定義されているVERBOSEに、-DVERBOSEを 指定してからmakeするとverboseモードになり、何も指定せずに makeするとopmは静かになります。 デフォルトではverboseモードになっています。

plotモードのインストール

結果を図示するプログラムを使いたいときには

  1. make plot

を行ってください。 ただし、PGPLOTが使える必要があります。PGPLOTは こちらで 入手することができます。同定を行って座標の変換式を求めたりするだけならば、 これらのプログラムは必要ありません。

make plotがうまく行えなかったあなた、ヘッダファイルが読込めなかったか、 PGPLOTライブラリとうまくリンクできなかったのかも知れません。 MakefileにあるPGDIRとCPGLIBSを適当に書き換えてください。

cfitsioモードのインストール

求めた座標の変換式をFITS画像のWCS(World Coordinate System)に 埋め込むことができます。本ソフト中では cfitsioのライブラリを利用して、ヘッダを編集します。 cfitsioが利用できる方は、MakefileのFITSDIR(fitsio.hが存在するディレクトリ)と FITSLIBS(ライブラリの呼び方)の設定を確認の上、

  1. make cfitsio
を行ってください。 具体的な使い方はこちらに説明があります。

各実行ファイルの説明

各実行ファイルの説明 (plotモード)

以下は、結果を図示するためのプログラムで、 PGPLOT が必要です。コンパイル方法については、こちら を参照してください。

各実行ファイルの説明 (cfitsioモード)

以下は、求めた変換式からWCSを埋め込むためのプログラムで、 cfitsio が必要です。コンパイル方法については、こちら を参照してください。なお、FITSへの埋め込みは行わずに、 必要なヘッダのキーワードと値だけを表示するプログラム (calc_wcs_header)もあります。

使い方の説明

コンパイルに成功したら、example1のディレクトリにいってみてください。 input1とinput2という２つのファイルがあります。 これは、2つの天体リストから「x y mag」を切り取ってきたものです。 input1の方はある画像のピクセル座標と機械等級、input2の方は カタログから取ってきた赤経赤緯とすでに補正された等級が並んでいます。 input1の座標がinput2の座標と合うように変換を行うので、 以下ではinput1を対象入力ファイル、input2を基準入力ファイルと呼びます。

まず、local_coordinate_RADecを使って、 赤道座標を局所平面座標に変換します。

../local_coordinate_RADec input2 input2b 17:46:10.6 -28:53:48.1 240

このデータの中心赤道座標と、視野サイズ(約480秒)の半分の値を与えました。 input2の座標を変換したinput2bという新しいファイルがつくられました。 これで、input1とinput2bを比較することができます。

次にopmを動かすときに必要なパラメータファイルをつくります。 ここではパラメータファイルの名前をparamとします。

../get_default_param param

このparamファイルにはデフォルトのパラメータが入っていますので、 必要な変更を行います。(適当なエディタを使っください。)

それでは、opmプログラムを動かしましょう。

../opm input1 input2b output param

計算はすぐに終わって、outputファイルが出力されます。 outputにはここに書いてあるような、 データが並びます。また、ファイルの最初の部分には、 使用したパラメータの値と、変換式が #で始まるコメント行として書かれています。 verboseモードならば、 モニターにもいろいろなログが 出力され、最後の部分に変換式が出ます。

plotモードのプログラムまでコンパイルしてあれば、 このoutputファイルを使っていろいろな図をプロットすることができます。

../plot_dxdy output 1
../plot_magdel output
../plot_map output 0 0 60

プロットの結果については、こちらをご覧下さい。

また、calc_wcs_headerを使えば、 求めた座標変換式に対応するようなWCSのキーワードと値をモニターに 表示することができます。ただし、この例で考えているように input1がピクセル座標で、 input2が局所球面座標を秒角単位で表している必要があります。

../calc_wcs_header output 17:46:10.6 -28:53:48.1

とすると、以下のような結果が表示されます。 また、cfitsioモードの put_wcs_headerを使えば、 これらのヘッダをFITS画像に直接書き込むことができます。 
WCSAXES  2
CTYPE1   'RA---TAN'
CTYPE2   'DEC--TAN'
CRPIX1   552.8094323
CRPIX2   489.5252822
CD1_1    -1.2666667e-04
CD1_2    1.1111111e-06
CD2_1    1.3888889e-06
CD2_2    1.2583333e-04
LONPOLE  180
CRVAL1   266.5441667
CRVAL2   -28.8966944
CUNIT1   'deg     '
CUNIT2   'deg     '

結果ファイルの書式

opmの結果得られるファイルは

id1 x1org y1org x1 y1 m1 id2 x2 y2 m2 flag

という書式です。それぞれの値の意味は以下の通り。

plotモードで描ける図

使い方の説明 でのデータを例にして、出力結果の図示を行いました。

cfitsioモードによるWCSの入力

赤道座標とピクセル座標との変換式がわかれば、そのFITS画像の WCS (World Coordinate System)がわかります。 そこで、put_wcs_headerを使えば、 opmの出力ファイルから直接FITSへのWCSの書き込みが行えます。 ただし、本ソフトは心射図法(gnomonic)によって局所球面座標から 射影平面座標への変換が行われていると仮定します。 また、input1(対象入力ファイル)はピクセル座標、input2(基準入力ファイル)は 局所球面座標が 秒角単位で書かれているようにしてください。局所球面座標は local_coordinate_RADecやlocal_coordinate_RADecを使えば、 赤経赤緯から秒単位で求めることができます。

たとえば、使い方の説明 は上記の条件を満たしています。対応する画像の名前が***.fitsだとしたら、

put_wcs_header ***.fits output 17:46:10.6 -28:53:48.1

とすれば、必要なヘッダが書き込まれます。なお、 calc_wcs_headerを使えば、 使い方の説明に書いた通り FITSIOによる書き込みは行わずに、 必要なキーワードと値をモニターに表示することができます。

使い方の説明 - 2. 合成カラー画像の作成

次の例として、3枚の画像について座標系を求めて、合成カラー画像を 作ってみましょう。そのために cfitsioモードのプログラム put_wcs_headerを利用します。 また、カラー画像を表示するためのソフトとして、 ds9を使います。 FITSIOをインストールしていない場合は、calc_wcs_header でWCSのパラメータを計算させてそれを何か別の方法でFITS画像に埋め込む こともできます。また、ds9は他のFITS画像表示ソフトでも代用できるでしょう。

example2のディレクトリに行くと、j.fits, h.fits, k.fitsという3枚の FITS画像が置いてあります。 それぞれJ(約1.2μm) H(約1.6μm) K(約2.2μm)という 3つのフィルターを使ってとった画像です。 そして、j.xym, h.xym, k.xymが写っている星のXY座標と等級、 catalog.xymが対応する領域のカタログ(2MASS)の局所平面座標(単位・秒)と Hバンド等級のファイルです。 catalog.xymの中心座標は、17:47:11.13 -26:49:10.4 (J2000)に対応しています。

以下のように、opmコマンドと put_wcs_headerコマンドを使うと、 各画像に必要なWCS情報が入ります。これをds9のRGBモードで表示すると カラー画像が見られます。(残念ながらあまりカラフルな星の無い領域ですが)

../opm j.xym catalog.xym j.match param
../opm h.xym catalog.xym h.match param
../opm k.xym catalog.xym k.match param
../put_wcs_header j.fits j.match 17:47:11.13 -26:49:10.4
../put_wcs_header h.fits h.match 17:47:11.13 -26:49:10.4
../put_wcs_header k.fits k.match 17:47:11.13 -26:49:10.4
ds9 -minmax -log -rgb -red k.fits -green h.fits -blue j.fits &

OPMパラメータの説明

opmプログラムでは、以下のパラメータを入力ファイルにいれて 動作させる。get_default_paramによって デフォルトのパラメータファイルをつくることができる。 カッコの中はデフォルトで設定される値。

自分のデータを使って同定するには

「x y mag」という書式の２つのファイルを準備してください。 x yは局所平面座標にしてください （ピクセル座標かだいたい共通する視野中心への接平面へ投影した座標）。 赤道座標になっている場合は、 local_coordinate_RADecやlocal_coordinate_degree を利用すれば、局所平面座標が得られます。 次にパラメータファイルを用意してください。 get_default_paramを使えば、デフォルトの 値を集めたパラメータファイルが作られますので、それを編集してください。 特に、MIRROR, MAGNIFY, XCENTER, YCENTER, HALFSIZE, TOLERANCEなどの パラメータに注意してください。

自分のデータをそのままの形で入力したいとき (上級編)

OPMプログラムの通常の入力は「x y mag」が各行に一つずつというものですが、 自分が持っているファイルをそのまま入力させたいとか、等級の誤差などの いろいろな情報も同じプログラムの中で処理させたいということがあります。 そのときは、original_data.hとoriginal_data.cを編集して、opm_original.cを 利用することを考えてください。 (C言語について、ポインタの扱いや分割コンパイルなどにあまり慣れていない人は、 通常のopmによる出力ファイルを適当に処理した方が容易でしょう。 system関数やシェルスクリプトを書けば、いくらでも処理を自動化できます。) 自分のファイルを入力する例を、example3のディレクトリで 試すことができます。とりあえず、OPMディレクトリでopm_originalを コンパイルしてから、example3でopm_originalを実行してみてください。

  1. make original
  2. cd example3
  3. ../opm_original input1 input2 output param

ここで与えられているinput1は通常の入力形式である「x y m」という 書式のファイルですが、input2はもっといろいろな情報を各行各天体に対して 含まれています。 

(input2の先頭５行)
17455437-2856402  -213.0  -172.2   8.542 0.027   6.889 0.046   6.193 0.026  AAA
17462835-2854488   233.1   -60.9   9.086 0.032   7.279 0.034   6.516 0.026  AAA
17461524-2850035    61.0   224.5  11.828 0.023   8.921 0.044   7.291 0.027  AAA
17460427-2857365   -83.0  -228.5  12.846 0.023   9.161 0.022   7.371 0.023  AAA
17460419-2856511   -84.1  -183.0   8.930 0.034   7.918 0.055   7.421 0.018  AAA
先頭から、ID(赤経赤緯を名前にしたもの)、局所平面座標でのx, y(秒角単位)、 Jバンド等級、Jバンド誤差、Hバンド等級、Hバンド誤差、Kバンド等級、Kバンド誤差、 JHKバンドそれぞれの精度を表すフラッグという書式です。このファイルをそのまま opmに入力することはできませんが、opm_originalを使うことで そのまま入力してしまい JHK３バンドの等級を出力ファイルに含めることができます。 (ただし、この例のoutputは通常の出力形式と異なるので、plotモードのプログラムに そのまま入力することはできません。)

このプログラムのメインはopm_original.cですが、original_data.[c,h]で 独自のファイル形式の入出力を加えている以外は、通常のメインプログラムopm.cと ほとんど同じです。original_data.[c,h]には、独自の書式のファイルを 入力したり(load_original_data)、そのデータをopmで操作する構造体(xym)に 変換したり(set_original_xym)、独自の形式で出力したり (output_matched_result_original)といった関数が記述されています。 細かい説明は省略しますが、独自のデータ形式へのポインタを保持するために、 xym構造体には、(void *)型のdata変数が定義されています。 xym構造体の中身はxym.c以外のソースから直接見ることはできませんが、 set_xym関数で入力し、get_xym_data関数でdata変数を取り出すことができます。 original_data.[c,h]を編集することによって、入出力の形式を 自分で好きなような形に変えられるはずです。opm_original.cとoriginal_data.[c,h] はそれぞれバックアップがありますので、自分のデータに あった入出力や処理を加えたいときには、これらのソースファイルを 編集してみて下さい。

トラブルシューティング

データの性質によっては、同定するのが難しい場合があります。 以下では、いくつか考えられる原因とその対処法を説明します。

OPMアルゴリズムについて簡単な説明

本同定ソフトは、Taburによって提案されたOptimistic Patter Matching (OPM) というアルゴリズムに基づいています。この方法の大まかな流れは以下の通りです。

  1. データの読込とソート 入力ファイルから「x y mag」を読み込んで、等級の明るい方から並べる。
  2. 三角形群の構築 等級の明るい方から、NSET個の星を使って様々な三角形をつくる。 そのぞれぞれについて、最長辺と最短辺の比(ratio)と、一番小さい角をはさむ二辺に よる内積(metric)を計算する。2つの入力ファイルがあるので、 三角形のグループが2つできる。
  3. 三角形マッチ : つくられた三角形に対して、互いのグループから１つずつ三角形を出し合って、 ratioとmetricが一致するペアを選びだす。
  4. Preliminary Verification : ratioとmetricの一致する三角形のペアが見つかったら、そのペアが 同じ３つの星からできていると仮定して、互いの座標系の変換式を計算してみる。 もし本当に同じ３つの星だったら、変換式は平行移動と回転だけで表せる １次式になっているはずなので、それを確認する。
  5. Final Verification : 前ステップで求めた変換式が適当な１次式だったら、 すべての星(あるいはNVERIFY個の明るい星同士)を使って、 その変換式に基づいた同定を行ってみる。 あらかじめ設定した割合(MATCHING)以上の星が ちゃんと両方の入力ファイルで同定できれば、その変換式が 正しいとして採用し、同定作業を終える。
座標を変換したあとに2つのファイルの天体を同定するときには、 以下のような方法でペアを探します。まず、対象入力ファイルも 基準入力ファイルも等級の明るい方から並び替えます（この並び替えは プログラムの中で行います）。 そして、対象入力ファイルの明るい方から同定天体を探していくのですが、 ある半径(TOLERANCE)に入るもっとも明るい天体を選んで、その天体は それ以降のリストから除きます。そうすると、お互いに近くにある 天体同士では明るい天体同士のペアから選ばれてリストからはずされていきます。 (この作業は２つのリストで明るさが同じとは限らないときに、予想外の カップリングを行う可能性があります。)

Taburの論文と本ソフトとで異なる点

本ソフトは、基本的にTaburの論文のアルゴリズムを採用していますが、 いくつか元の論文の記述とは異なる方法で計算している部分があります。

著作権や利用の条件について

OPM: a pattern matching software
Copyright (c) 2007-2008 Noriyuki Matsunaga

本ソフトウェアは、学術的な目的において、 個人・団体によらずに使用・改良・再配布を許可します。 ただし、再配布を行う場合は、オリジナルの内容が それとわかる状態で一緒にダウンロードされるようにして、 本利用条件を変更せずにReadMeなどへ提示してください。 商用目的のためには、原則的にすべての内容について 利用を禁止します。もしそのような希望があれば、松永までご相談下さい。 著作権は放棄しておらず、松永典之個人に帰属します。

本ソフトは、Tabur (2007)の記述に基づき、 松永個人が作成したものです。最低限の動作確認は しているつもりですが、自己責任で使用してください。 本ソフトのコンパイル・使用などすべての動作で起きうる不利益に関する 一切の責任を放棄します。 また、Taburが記述した通りに動作することを 保証するものではありません。

また、本ソフトに含まれるコードのうち、thdrad.ctmat.cwcs.cは 田辺俊彦氏が作成したものを再利用しました。したがって、 それらに対する著作権は田辺氏にあります。田辺氏のご好意により、 本ソフトで利用し、他のファイルと同様に公開できることを この場を借りて感謝いたします。なお、それらのプログラムが、 本ソフト内で正しく利用され、本来の機能を発揮しているか どうかの判断は、松永の責任において行いました。
thdrad.c, tmat.c, wcs.c: Copyright (c) 2007 Toshihiko Tanabe

ソフトウェアの性能に関わるような改良を 加えられた場合は、ぜひとも 松永 までご連絡ください。

引用・謝辞の入れ方

論文などで引用、あるいは謝辞に入れてくださる場合は、 以下の例を参考にしてください。

Copyright (c) 2007-2008 Noriyuki Matsunaga