SLLIB,SFITSIO は宇宙研の赤外の皆様,宇宙研のデータセンターの皆様の 貴重な協力と支援のもので開発されてきたソフトウェアです.
特に宇宙研のデータセンターでは,これらのソフトウェアの価値に 多大なるご理解をいただき, 優秀な人材と多額な資金を惜しみなく SLLIB と SFITSIO の開発 (試験・マニュアル作成・翻訳) に投入してくださいました.
感謝を表明したいと思います.
データセンター・解析プロジェクトなどでの業務や, 真面目にデータ解析を行ないたい天文学者のための, 「FITS I/O 環境」と「スクリプト言語のようなプログラミング言語環境」を提供する, C/C++ のライブラリ・セット.
メモリ管理,配列の選択・コピーなどの 「相対的に科学と距離のある部分」に関する数々の無駄な労力を除去し, プログラマが 科学の本質部分のコーディングに集中でき, 「正しさ」と「透明さ」を保ちつつ, 高いパフォーマンスを狙える開発環境 の整備を目的としている.
printf("telescop = %s\n",
fits.image("Primary").header("TELESCOP").svalue());
C++の言語仕様は優れているのだが…
スクリプト言語にも長所と短所があり, C/C++ に対する需要は無くなりそうにない.
次の2つを作れば良さそうである:
Makefile
の作成は不要で,
「s++ my_program.cc
」
でコンパイルは完了.
【例】 arr1 = arr0.sectionf("0:99,*");
【例】 arr1 = arr0.sectionf("%d:%d,*", begin, end);
.gz
,.bz2
)
を直接読み書き可能.【例】 fits.read_stream("http://foo.com/myimage.fits.gz");
【例】 fits.read_stream("myimage.fits.gz[1:100,*]");
CONTINUE
,CHECKSUM
,DATASUM
キーワード等,広く使われている非標準FITS拡張に対応.
【例】 f_in.open("r", "http://www.xxx.jp/datafile.dat.gz");
【例】 tstring my_url = "http://darts.isas.jaxa.jp/foo/";
my_url.regreplace("([a-z]+://)([^/]+)(.*)", "\\2");
sli/complex.h
).宇宙研ではすでにアーカイブの業務や衛星プロジェクトでの基幹システムにて, 採用実績を有する.
s++
」コマンドで,
テンプレートコードを作る.$ s++ fits_read_write.cc -lsfitsio
Do you want to create a template source file? [y/n] y A template code is written: fits_read_write.cc
#include <sli/stdstreamio.h> #include <sli/fitscc.h> using namespace sli; /* * main function */ int main( int argc, char *argv[] ) { stdstreamio sio; fitscc fits; sio.printf("Hello World\n"); return 0; }
s++
」コマンドで
コンパイル(「-o」等はいらない).$ s++ fits_read_write.cc -lsfitsio g++ -I/usr/local/include -L/usr/local/lib -Wall -O2 -o fits_read_write fits_read_write.cc -lsfitsio -lsllib -lreadline -lcurses -lbz2 -lz
$ s++ fits_read_write.cc -lsfitsio / fitsfile.fits
#include <sli/stdstreamio.h> #include <sli/digeststreamio.h> 他
stdstreamio sio; /* 標準出力用のオブジェクト */ sio.printf("標準出力\n"); sio.eprintf("標準エラー出力\n");従来の
printf(…);
, fprintf(stderr,…);
を使っても OK.
digeststreamio f_in; f_in.open("r", "http://foo.ac.jp/file.txt.gz");「digeststreamio」は様々な圧縮・非圧縮ストリームに対応する最強のクラス
for ( i=0 ; i < N ; i++ ) { f_in.openf("r", "ftp://foo.ac.jp/file_%d.txt.gz", i); : f_in.close(); }
const char *line; while ( ( line=f_in.getline() ) != NULL ) { printf("%s",line); }⇒ デモ
#include <sli/tstring.h>
stdstreamio sio; tstring my_str; /* 文字列オブジェクトを作る */ my_str = "Hello World"; /* "Hello World" を my_str に代入 */ sio.printf("my_str = '%s'\n", my_str.cstr());
my_str = "Hello World"; my_str[13] = '!'; my_str[14] = '!'; sio.printf("my_str = '%s'\n", my_str.cstr());《実行結果》
my_str = 'Hello World !!'
my_str.printf("%s %s", "Hello", "World"); /* 上書き */ my_str.appendf(" %s", "!!"); /* 追加 */
my_str = " Hello World "; sio.printf("before: '%s'\n", my_str.cstr()); my_str.trim(); sio.printf("after trim(): '%s'\n", my_str.cstr()); my_str.tolower(); sio.printf("after tolower(): '%s'\n", my_str.cstr());《実行結果》
before: ' Hello World ' after trim(): 'Hello World' after tolower(): 'hello world'
sed -e …
」と同じ事をやってみる)
stdstreamio sio; tstring my_url = "http://darts.isas.jaxa.jp/foo/"; my_url.regreplace("([a-z]+://)([^/]+)(.*)", "\\2", false); sio.printf("hostname = %s\n", my_url.cstr());《実行結果》
hostname = darts.isas.jaxa.jp
stdstreamio sio; /* 0123456789012345678901234567890 */ tstring str = "hypot( creal(v), cimag(v) )"; size_t span; ssize_t pos = str.find_quoted("()", '\\', &span); sio.printf("pos = %zd span = %zu\n", pos, span);《実行結果》
pos = 5 span = 22
#include <sli/tarray_tstring.h>
tarray_tstring my_arr; my_arr[2] = "fuji"; my_arr[3] = "hayabusa"; my_arr.dprint();《実行結果》
sli::tarray_tstring[obj=0xbfffee08] = {"", "", "fuji", "hayabusa"}「
.dprint()
」はオブジェクトの内容を表示するもので,
デバッグに使う.Python の「print
」,
PHP の「print_r()
」,IDL の「help,
」のような
もの.
argv[]
を代入
int main( int argc, char *argv[] ) { tarray_tstring args = argv;
const char *line = " Z-80,, 8086 , 6800"; tarray_tstring my_arr; my_arr.split(line, ",", true); my_arr.dprint();《実行結果》
sli::tarray_tstring[obj=0xbfffee28] = {" Z-80", "", " 8086 ", " 6800"}「
.split()
」の3つめの引数が false の場合,
文字列長ゼロの要素が作られない.
2つめの引数は「区切り文字」の文字集合であり,
" \t"
のように指定できる.
my_arr.trim(); my_arr.dprint();《実行結果》
sli::tarray_tstring[obj=0xbfffee28] = {"Z-80", "", "8086", "6800"}
stdstreamio sio; tarray_tstring my_elms; tstring my_str = "OS = linux"; my_elms.regassign(my_str, "([^ ]+)([ ]*=[ ]*)([^ ]+)"); if ( my_elms.length() == 4 ) { sio.printf("keyword=[%s] value=[%s]\n", my_elms.cstr(1), my_elms.cstr(3)); }《実行結果》
keyword=[OS] value=[linux]
const char *line = "winnt( ) 'program files'"; tarray_tstring my_arr; my_arr.split(line, " ", false, "'()", '\\', false); my_arr.dprint();《実行結果》
sli::tarray_tstring[obj=0xbfffee28] = {"winnt( )", "'program files'"}
内部では単純な1次元バッファを持つ. メンバ関数の引数は 0-indexed である.
型も配列長も可変だが, 型に依存したメンバ関数が使えない.
通常はこちらを使う. 型が固定である事以外のデメリットはない. 異なる型どおしでも,演算が可能.
#include <sli/mdarray.h> /* mdarrayクラスとその継承クラス */ #include <sli/mdarray_statistics.h> /* 統計用の関数 */ #include <sli/mdarray_math.h> /* 配列に対する数学関数 */
/* オブジェクトを作成.*/ /* 引数は assign() 等での自動リサイズ on/off を設定 */ mdarray_float arr0(false); /* 配列の大きさを設定し,値を代入 */ arr0.resize_2d(8,4); arr0 = 500.0; arr0 *= 2; arr0.dprint();《実行結果》
sli::mdarray[obj=0xbffff4e0, sz_type=-4, dim=(8,4)] = { { 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000 }, { 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000 }, { 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000 }, { 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000 } }
size_t i, j; for ( i=0 ; i < arr0.row_length() ; i++ ) { /* Y */ for ( j=0 ; j < arr0.col_length() ; j++ ) { /* X */ arr0(j,i) = 100 + 10*i + j; } } arr0.dprint();
arr0(x,y)
の引数は 0-indexed である.
float *const *arr0_ptr = arr0.array_ptr_2d(true); size_t i, j; for ( i=0 ; i < arr0.row_length() ; i++ ) { /* Y */ for ( j=0 ; j < arr0.col_length() ; j++ ) { /* X */ arr0_ptr[i][j] = 100 + 10*i + j; } } arr0.dprint();《実行結果》
sli::mdarray[obj=0xbffff4e0, sz_type=-4, dim=(8,4)] = { { 100, 101, 102, 103, 104, 105, 106, 107 }, { 110, 111, 112, 113, 114, 115, 116, 117 }, { 120, 121, 122, 123, 124, 125, 126, 127 }, { 130, 131, 132, 133, 134, 135, 136, 137 } }「
.array_ptr_2d()
」はオブジェクト内部で作成される
ポインタ配列(アドレステーブル)を返す.
このポインタ配列はリサイズ等が行なわれると,自動的に更新される..array_ptr_3d()
」もある.
mdarray_double arr1(false); arr1 = arr0.sectionf("1:4, *"); arr1.dprint();《実行結果》
sli::mdarray[obj=0xbffff2c0, sz_type=-8, dim=(4,4)] = { { 101, 102, 103, 104 }, { 111, 112, 113, 114 }, { 121, 122, 123, 124 }, { 131, 132, 133, 134 } }「
arr1 = arr0.sectionf("1:4, *");」
の「=
」
では,shallow copy が使われる."1:4,*"
または "(1:4,*)"
の場合,0-indexed である.ただし,"[…]"
の場合は
1-indexed であり,
今の場合は "[2:5,*]"
とも書ける.
arr0.pastef(arr1, "*, 2:3"); arr0.dprint();《実行結果》
sli::mdarray[obj=0xbffff4e0, sz_type=-4, dim=(8,4)] = { { 100, 101, 102, 103, 104, 105, 106, 107 }, { 110, 111, 112, 113, 114, 115, 116, 117 }, { 101, 102, 103, 104, 124, 125, 126, 127 }, { 111, 112, 113, 114, 134, 135, 136, 137 } }
arr0.dividef(arr1, "*, 2:3"); arr0.dprint();《実行結果》
sli::mdarray[obj=0xbffff4e0, sz_type=-4, dim=(8,4)] = { { 100, 101, 102, 103, 104, 105, 106, 107 }, { 110, 111, 112, 113, 114, 115, 116, 117 }, { 1, 1, 1, 1, 124, 125, 126, 127 }, { 1, 1, 1, 1, 134, 135, 136, 137 } }
arr0 = pow(arr0, 2.0); arr0.dprint();《実行結果》
sli::mdarray[obj=0xbffff4e0, sz_type=-4, dim=(8,4)] = { { 10000, 10201, 10404, 10609, 10816, 11025, 11236, 11449 }, { 12100, 12321, 12544, 12769, 12996, 13225, 13456, 13689 }, { 1, 1, 1, 1, 15376, 15625, 15876, 16129 }, { 1, 1, 1, 1, 17956, 18225, 18496, 18769 } }
/* get mean, variance, skewness, kurtosis */ mdarray_double moment = md_moment(arr0, false, NULL, NULL); moment.dprint();《実行結果》
sli::mdarray[obj=0xbfffee80, sz_type=-8, dim=(4)] = { 10165.625, 41496969.79, -0.6201298573, -1.038419883 }
arr1 = md_median_y(arr0); arr1.dprint();《実行結果》
sli::mdarray[obj=0xbffff2c0, sz_type=-8, dim=(8,1)] = { { 5000.5, 5101, 5202.5, 5305, 14186, 14425, 14666, 14909 } }
/* 3次元アレイを用意 */ arr0.resize_3d(1024, 512, 10); /* 1枚ずつ z の位置を変えて 3次元アレイにペースト */ mdarray_float arr_tmp; arr_tmp.resize_2d(1024, 512); for ( i=0 ; i < 10 ; i++ ) { : /* arr_tmp に代入する処理 */ : arr0.pastef(arr_tmp, "*,*,%d", i); } /* z 方向に median コンバイン */ arr1 = md_median_small_z(arr0);コンバインする画像が多い場合は,
md_median_small_z()
のかわりに,
md_median_z()
を使う事ができる.
前者は z-x 面単位で一時バッファを使うので,
高速だがメモリ使用量が多い.
後者は z 方向の1本単位で一時バッファを使うので,
省メモリだが,z が小さい場合にパフォーマンスが出ない.
画像,テーブル,ヘッダ,HDUとも, メンバ関数の引数は常に 0-indexed である (ファイル読み取り時の範囲指定時のみ,1-indexed でも記述可能).
#include <sli/fitscc.h> #include <sli/fits_image_statistics.h>
int main( int argc, char *argv[] ) { fitscc fits; fits.read_stream(argv[1]); fits.write_stream(argv[2]); }《使用例》
$ ./fits_read_and_write org_image.fits.gz'[1:100,*]' trim_image.fits
[1:100,*]
の場合は 1-indexed である.
0-indexed で書きたい場合は,
(0:99,*)
と書く.printf("N_HDU = %ld\n", fits.length());
/* read */ double gain = fits.image(0L).header("GAIN").dvalue(); /* write */ fits.image(0L).header("CRVAL1").assign(1.234);
.dvalue()
… double型の値を取得.lvalue()
… long型の値を取得.svalue()
… 文字列(const char *
)型の値を取得
image(0L)
」は,
「image("Primary")
」でも良い.
fits_image &prim = fits.image(0L);「
fits.
」から書くと長いので,
短縮表記できるようにしている.
if ( 0 < prim.header_value_length("EQUINOX") ) { /* OK */ }キーワードが存在しない場合は,返り値が負となり, 値が存在しない場合は,返り値が 0 となる.
int r_type = prim.header("EQUINOX").type(); if ( r_type == FITS::DOUBLE_T ) { /* 実数 */ } else if ( r_type == FITS::LONGLONG_T ) { /* 整数 */ } else if ( r_type == FITS::DOUBLECOMPLEX_T ) { /* 複素数 */ } else if ( r_type == FITS::BOOL_T ) { /* 論理値 */ } else if ( r_type == FITS::STRING_T ) { /* 文字列 */ } else { /* 予期せぬエラー */ }
double v = prim.dvalue(0,99); prim.assign(v, 1,99);引数は 0-indexed の x, y の順で,速度は遅い. fixpix のようなものを作る時はこれで十分.
prim.convert_type(FITS::FLOAT_T);
mdarray_float im_arr(false); prim.data_array().swap(im_arr);速度や汎用性のためには,mdarray_* のオブジェクトにデータを移し, そちらのメンバ関数・関数を使うと良い.
im_arr *= gain; mdarray_double moment = md_moment(im_arr.sectionf("2:21, *")); moment.dprint();
mdarray_double moment = fitsim_moment(prim.sectionf("2:21, *"); moment.dprint();
「FITSテンプレート」とは,新規 FITS の元となる, FITSヘッダのようなテキストファイルである. SFITSIO は,曖昧な FITS テンプレートから, 正しい FITS を生成してくれる.
# # Test for Primary only HDU # NAXIS2 = 512 NAXIS1 = 512 BITPIX = 32 COMMENT ---------------------------------------------------------------- EQUINOX = 2000.0 CTYPE1 = RA---TAN CTYPE2 = DEC--TAN CRPIX1 = 0.0 CRPIX2 = 0.0 CRVAL1 = 0.0 CRVAL2 = 0.0 CDELT1 = 0.1 CDELT2 = 0.1 PC1_1 = 1.0 PC1_2 = 0.0 PC2_1 = 0.0 PC2_2 = 1.0
fitscc fits; fits.read_template("template.txt"); fits.write_stream("my_image.fits.gz");
SIMPLE = T / conformity to FITS standard BITPIX = 32 / number of bits per data pixel NAXIS = 2 / number of data axes NAXIS1 = 512 / length of data axis 1 NAXIS2 = 512 / length of data axis 2 EXTEND = F / possibility of presence of extensions COMMENT ---------------------------------------------------------------- EQUINOX = 2000.0 / equinox of celestial coordinate system CTYPE1 = 'RA---TAN' / type of celestial system and projection system CTYPE2 = 'DEC--TAN' / type of celestial system and projection system CRPIX1 = 0.0 / pixel coordinate at reference point CRPIX2 = 0.0 / pixel coordinate at reference point CRVAL1 = 0.0 / world coordinate at reference point CRVAL2 = 0.0 / world coordinate at reference point CDELT1 = 0.1 / world coordinate increment at reference point CDELT2 = 0.1 / world coordinate increment at reference point PC1_1 = 1.0 / matrix of rotation (1,1) PC1_2 = 0.0 / matrix of rotation (1,2) PC2_1 = 0.0 / matrix of rotation (2,1) PC2_2 = 1.0 / matrix of rotation (2,2)
データセンターなどでは,大量の FITS ファイルを管理するために, FITS ヘッダを高速に読み取る必要がある.
FITS ヘッダをマネージするための fits_header クラスが持つ
read_stream()
を使えば,
Data Unit を読まずにファイルへのアクセスを終了する事が可能.
もちろん,圧縮ファイルでも高速アクセスが可能.
#include <sli/fitscc.h> #include <sli/digeststreamio.h> using namespace sli;
.read_stream()
でヘッダ部分のみ読み取る.
int main() { digeststreamio f_in; fits_header hdr; /* ヘッダ部分のみ読む */ f_in.open("r", "my_image.fits.gz"); hdr.read_stream(f_in); /* ヘッダの内容を表示 */ printf("gain = %g\n", hdr.at("GAIN").dvalue()); printf("[all string]\n"); printf("%s\n", hdr.formatted_string());この時点でファイルポインタは次のData Unitの先頭にあるので,
f_in.read(...)
などでData Unitにアクセス可能.
hdr.skip_data_stream(f_in);
hdr.read_stream(f_in);
f_in.close();