次のコードは，プログラム側で sin() と cos() の関数の値を計算し， gnuplot に FIFO でデータを与えて描画させるためのコードです． 極めて平凡なコードである事がおわかりいただけると思います．

#include <math.h>
#include <sli/stdstreamio.h>
#include <sli/pipestreamio.h>
using namespace sli;

int main( int argc, char *argv[] )
{
    int ret_status = -1;
    pipestreamio p_out;
    stdstreamio sio, f_out;
    int i;

    /* gnuplot を起動 */
    if ( p_out.open("w", "gnuplot") < 0 ) {
        sio.eprintf("[ERROR] p_out.open() failed\n");
        goto quit;
    }

    /* gnuplot へコマンドを送信 */
    p_out.printf("plot \"fifo0\" with lines, \"fifo1\" with lines\n");
    p_out.flush();

    /* fifo0 経由で gnuplot へデータを送信 */

    if ( f_out.open("w","fifo0") < 0 ) {
        sio.eprintf("[ERROR] f_out.open() failed\n");
        goto quit;
    }

    for ( i=0 ; i < 1000 ; i++ ) {
        double x = 4 * M_PI * i / 1000.0;
        f_out.printf("%.15g %.15g\n", x, sin(x));
    }

    f_out.close();

    /* fifo1 経由で gnuplot へデータを送信 */

    if ( f_out.open("w","fifo1") < 0 ) {
        sio.eprintf("[ERROR] f_out.open() failed\n");
        goto quit;
    }

    for ( i=0 ; i < 1000 ; i++ ) {
        double x = 4 * M_PI * i / 1000.0;
        f_out.printf("%.15g %.15g\n", x, cos(x));
    }

    f_out.close();

    /* ターミナルから入力を待つ */
    sio.getchr();
    
    /* gnuplot への接続をクローズ */
    p_out.close();

    ret_status = 0;
 quit:
    return ret_status;
}

注意点は，FIFO に書き込む順番は， gnuplot の plot コマンドに与えた FIFO の順に一致させる事です．

○実習10○

1. 上記のコードの動作確認を行なってください．
2. プロットする関数を変えてみてください．

参考:
SLLIB ダイジェスト版 リファレンス 「ストリーム」

文字列

SLLIB では文字列を扱うのに 「tstring」 というクラスを使います． tstringクラス はスクリプト言語のように文字列を扱えるようにするもので， オブジェクトの内部に自動管理されるヒープバッファを持ち，そこに 文字列が記録されます． バッファの大きさの変更，バッファの開放は自動的に行なわれ， 高い安全性を持ちます．

文字列の扱いの基本

早速コードを見ていきましょう．

#include <sli/stdstreamio.h>
#include <sli/tstring.h>
using namespace sli;

int main( int argc, char *argv[] )
{
    stdstreamio sio;
    tstring my_str;                        /* 文字列オブジェクトを作る */
    my_str = "Hello World";                /* "Hello World" を my_str に代入 */
    sio.printf("my_str = '%s'\n", my_str.cstr());

    return 0;
}

C言語で，

    const char *my_str = "Hello World";

とすると，文字列定数 "Hello World" の先頭アドレスが my_str に記録されましたが， tstring の場合は，文字列定数 "Hello World" の文字列が オブジェクト内部のバッファにコピーされます． そして，printf() の %s などで 文字列の先頭アドレスが必要な場合は， .cstr() で内部バッファのアドレスを取得します．

これだけのコードだと 一見面倒なように見えますが，文字列の変更のような事をしようとすると， C言語の時よりずっと簡単かつ安全に目的の処理を行なう事ができます．

○実習11○

1. 上記のコードについて，動作確認を行なってください．
2. 上記コードの return 文の前に，下記のコードを挿入し， その結果を確認してください．

       my_str[13] = '!';
       my_str[14] = '!';
       sio.printf("my_str = '%s'\n", my_str.cstr());
   

実習11 では， [i] でアクセスした場合も， 内部バッファの大きさは自動的に調整される事が理解できたと思います．

文字列の編集

tstring クラスは，バッファの自動化だけでなく， 様々な文字列編集のためのメンバ関数が揃っている事も特徴の1つです． 文字列の追加，文字列の左右の空白の除去，小文字への変換を行なってみましょう．

#include <sli/stdstreamio.h>
#include <sli/tstring.h>
using namespace sli;

int main( int argc, char *argv[] )
{
    stdstreamio sio;
    tstring my_str;                              /* 文字列オブジェクトを作る */

    my_str = "  Hello ";                         /* "  Hello " を代入 */
    my_str.append("WORLD ");                     /* "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());

    return 0;
}

○実習12○

上記のコードの動作確認を行なってください．

文字列の検索，拡張正規表現

単純な文字列比較，シェル風のマッチ，POSIX 拡張正規表現を使って， オブジェクト内文字列の検索，置換が可能です．

次の例は，URL 文字列からホスト名を取り出すのに， 正規表現を使ったものです． regreplace() メンバ関数は，拡張正規表現(第1引数)によるマッチを行ない，マッチした部分を 任意の文字列(第2引数)に置き換えるメンバ関数です． 次のように，sed コマンドと同様，「\\1」「\\2」等で 後方参照が使えます．

    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());

ストリームからのテキストファイルの読み出し

ストリームのところで登場した getline() メンバ関数と tstring クラスとを組み合わせると， 次のようにスクリプト言語に近い形で， テキストファイルの読み出しを行なう事ができます．

#include <sli/stdstreamio.h>
#include <sli/digeststreamio.h>
#include <sli/tstring.h>
using namespace sli;

/*
 * An example code for reading a text file from network
 */

int main( int argc, char *argv[] )
{
    int ret_status = -1;
    stdstreamio sio;                            /* stdin, stdout and stderr */
    digeststreamio f_in;                        /* for local file */

    if ( 1 < argc ) {

        const char *filename = argv[1];
        tstring line;                           /* buffer for a line */
        
        /* ストリームをオープン */
        if ( f_in.open("r", filename) < 0 ) {
            sio.eprintf("[ERROR] cannot open: %s\n",filename);
            goto quit;
        }

        /* テキストファイルを 1 行ずつ読む */
        while ( (line=f_in.getline()) != NULL ) {
            line.chomp();                               /* erase "\n" */
            sio.printf("%s\n",line.cstr());
        }

        /* ストリームをクローズ */
        f_in.close();

    }

    ret_status = 0;
 quit:
    return ret_status;
}

chomp() は，文字列の右端の改行文字を除去するためのメンバ関数です． 上記のコードの while() の部分は，SLLIB で良く使う コードのパターンですので，覚えておくと便利です．

参考:
SLLIB ダイジェスト版 リファレンス 「文字列」

文字列配列

SLLIB で文字列配列を扱うには， tarray_tstringクラス を使います． このクラスも tstring クラスと同様， ほぼ完全なメモリの自動管理により，易しく安全に文字列配列を扱う事ができま す．C言語のポインタ配列との親和性も考慮した設計と API により， 様々な場面で応用できるのも特徴の1つです．

文字配列の扱いの基本

tstring では， あらかじめバッファを確保していなくても [i] で いきなり代入できました．これは，tarray_tstring でも同じです． 例を示します．

#include <sli/stdstreamio.h>
#include <sli/tarray_tstring.h>
using namespace sli;

int main( int argc, char *argv[] )
{
    tarray_tstring my_arr;

    my_arr[2] = "fuji";
    my_arr[3] = "hayabusa";

    my_arr.dprint();

    return 0;
}

このコードの最後の dprint() は， オブジェクトの内容を標準エラー出力に出力するものです． 次のような表示を得ます．

sli::tarray_tstring[obj=0xbfffee08] = {"", "", "fuji", "hayabusa"}

これにより，いちいち for ループを書かなくても， 内容の確認を手早く行なう事ができます． もちろん，for ループで1つ1つ内容を表示する事もできて， その場合は，次のように書きます．

    stdstreamio sio;
    size_t i;
    for ( i=0 ; i < my_arr.length() ; i++ ) {
        sio.printf("my_arr[%zu] = '%s'\n", i, my_arr[i].cstr());
    }

length() で配列の個数， [i] で tstring オブジェクトの参照を取得します． 注意していただきたいのは，[i] の後についている .cstr() は tstring クラスのメンバ関数 であるという事です．これはつまり， tarray_tstring オブジェクトは， その内部に配列要素の個数分のtstringオブジェクトを持っていて， [i] で該当する内部の tstring オブジェクトを 引き出している，という事になります．

このように，オブジェクトの内部で，別のクラスの複数のオブジェクトを管理している， という手法はオブジェクト指向の基本的なパターンであり， 今後もずっと出てきますので，注意してみていてください． 内部のオブジェクトを引き出すメンバ関数の形としては， [i] や .at(i) が多いですが， より具体的な名前である事もあります．

○実習13○

1. 上記のコードの動作確認を行なってください．
2. 以下は csv 形式の文字列を分割するためのコードの一部です．

       const char *line = " Z-80,,  8086 , 6800";
       tarray_tstring my_arr;
       my_arr.split(line, ",", true);
       my_arr.dprint();
   

   動作確認を行なってください．
3. my_arr.split(line, " ", false); のように書くと， 空白区切りに対応します．動作確認を行なってください．
4. 2. で作ったコードに，次のコードを追加するとどのようになるか， 確認してください．

       my_arr.trim();
       my_arr.dprint();
   

char *argv[] の情報を格納

文字列のポインタ配列の情報を，まるごと配列にぶちこむのも = 記号 一発です．

#include <sli/tarray_tstring.h>
using namespace sli;

int main( int argc, char *argv[] )
{
    tarray_tstring args;

    args = argv;
    args.dprint();

    return 0;
}

ここから正規表現などを使って引数をチェックしたり， 値の変換をしたり，チェック済みの引数を削除したり， 様々な応用が安全かつ簡単に行なえます．

正規表現マッチの結果(後方参照の情報を含む)を格納

正規表現の (…) により， 文字列をいくつかのパーツへ分割したい事があります． このような場合は， regassign() メンバ関数を使うと， マッチした部分の情報を tarray_tstring の文字列配列として取得する事ができます． 例を示します．

    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[1].cstr(), my_elms[3].cstr());
    }

これを実行すると，次の結果を得ます．

keyword=[OS] value=[linux]

このように正規表現の (…) が3つのパートを指定し，それが マッチした場合は，文字列配列の長さが 4 になります． 配列の各内容については， my_elms[0] にマッチした全体の部分， my_elms[1] に1つめの (…) 部分， my_elms[2] に2つめの (…) 部分，のように セットされます．

高度な文字列分割 (括弧やクォーテーションを認識)

次のような高度な文字列分割も行なえます．空白区切りの文字列を 分割しますが， 括弧内の空白と，クォーテーション内の空白は区切り文字とはみなしません．

    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'"}

◎課題3◎

1. これまで出てきたコードを参考に， このテキストファイル を読み込み， 各行の処理において文字列配列の要素に値を取り込んでください． tstring クラスの chomp()， tarray_tstring クラスの split() を使ってみてください．
2. ra[deg], dec[deg] から 単位ベクトル cx, cy, cz を計算して 結果を表示してください．文字列から実数への変換は tstring クラスの atof() を使ってみてください．

文字列連想配列

ここまでくると，もうほとんどスクリプト言語と変わらなくなってきます． 文字列連想配列のためのクラスは asarray_tstring で， さきほどの文字列配列と同じように使えます．

#include <sli/stdstreamio.h>
#include <sli/asarray_tstring.h>
using namespace sli;

int main( int argc, char *argv[] )
{
    stdstreamio sio;
    asarray_tstring my_arr;                   /* 文字列の連想配列 */
    size_t i;

    /* 代入 */
    my_arr["JR-EAST"] = "SUICA";
    my_arr["JR-CENTRAL"] = "TOICA";
    my_arr["JR-WEST"] = "ICOCA";

    /* 表示 (もちろん .dprint() を使っても良い) */
    for ( i=0 ; i < my_arr.length() ; i++ ) {
        const char *key = my_arr.key(i);
        sio.printf("%s: [%s]\n", key, my_arr[key].cstr());
    }

    return 0;
}

文字列配列の場合と同様，内部で複数の tstring のオブジェクトを 管理しており，my_arr[キー文字列] の形で 該当する tstring のオブジェクトを引き出します．簡単ですね．

○実習14○

1. 上記のコードの動作確認を行なってください．

多次元配列

SLLIB で提供される『3点セット』--- 「1.ストリーム」 「2.文字列」 「3.多次元配列」 のうち， 多次元配列は，科学研究においては特に高い需要があります． にもかかわらず，C++ の世界には多次元配列を扱うのに適当なクラスは あまり見かけません．したがって， SLLIB の最大の強みは，この多次元配列にあり，と言えるでしょう．

SLLIB で多次元配列を扱うためのクラスは mdarray とその継承クラス であり， n次元の配列を最新のスクリプト言語と同様の感覚で扱う事ができます． メモリ管理の自動化はもちろん， mdarray オブジェクトの内部は常に単純な1次元バッファとなっている上， 内部バッファへの直接アクセスを前提とした安全設計により， 従来の C言語のライブラリとのリンクも容易です． 当然，配列の要素の数え方は他のクラスと同様 0-indexed です．

クラスの選び方

多次元配列のクラスの使い方はストリームの場合と似ており， 「基底クラス」「継承クラス」があります． 基底クラスは「mdarray」で， 継承クラスは「mdarray_double」「mdarray_float」のように 配列で扱う型の名称がつきます．また， 基底クラスは「抽象」ではないので，プログラマが使う事ができます．

通常，プログラマは次に示す継承クラスを使います． 型が固定である事以外のデメリットはなく，異なる型どうしでも， 四則演算を含む様々な演算が可能です．

配列データを「継承クラス」「基底クラス」の間で移動させる事ができるので， 必要になった時点でクラスを切り替える事も容易です．

多次元配列の基本

mdarray とその継承クラスには，2つの動作モード 「自動リサイズモード」「手動リサイズモード」があります． 前者はこれまで見てきた文字列配列のように，[i] のような 要素へのアクセスで自動的にバッファを再確保するモードです． 後者は std::vector のようにバッファの自動確保をせず resize() メンバ関数など を使ってプログラマがバッファの大きさを設定するモードです． どちらも，スコープを抜けた時のメモリの自動開放と， [i] のようなアクセスでの バッファ境界のチェックが入るなどの安全性は文字列配列などの場合と同じです．

観測データの場合は，2次元以上のデータが多く，バッファの大きさは 頻繁にリサイズしないと思いますので， 基本的には「手動リサイズモード」を使う方向で良いと思います．

次に，基本的なコードの例を示します．

#include <sli/mdarray.h>             /* mdarrayクラスとその継承クラス */
#include <sli/mdarray_statistics.h>  /* 統計用の関数 */
#include <sli/mdarray_math.h>        /* 配列に対する数学関数 */
using namespace sli;

int main( int argc, char *argv[] )
{
    /* オブジェクトを作成 */
    /* 引数は assign() 等での自動リサイズ on/off を設定 */
    mdarray_float arr0(false);

    /* 配列の大きさを設定し，値を代入 */
    arr0.resize_2d(8,4);
    arr0 = 500.0;
    arr0 *= 2;
    arr0.dprint();

    return 0;
}

まず，ヘッダファイルは sli/mdarray.h を include すれば すべての継承クラスが使用できるようになります．

オブジェクト arr0 を作るところで， 引数を与えていますが，これが動作モードを決めます．この例では， 手動リサイズモードを選択しています．その後， resize_2d() で 8 × 4 の2次元配列とし， 全配列要素に対してスカラー値 500.0 を代入，スカラー値 2 を掛けています． これを実行すると次のような結果を得ます．

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 }
}

○実習15○

上記のコードの動作確認を行なってください．

一要素ずつの代入

3つの方法を示します． まずは，安全な方法です．

    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) のようにアクセスします． この方法では常にバッファ境界のチェックが入り， 範囲外に書き込んでも単に無視されます (範囲外を読み込んだ場合は NAN が返ります)． 次に示す方法に比べて性能は劣りますが， スクリプト言語に比べれば十分高速です． 第3部では，この丸括弧で要素にアクセスする方法を使っています．

次は高速な方法です．

    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();

arr0.array_ptr_2d(true) とすると，2次元データにアクセスするためのポインタ配列を オブジェクト内部で生成してくれるので，それを使います． こうする事で， arr0_ptr[y][x] の形で各要素に高速にアクセスする事ができます． このポインタ配列はリサイズ等が行なわれると，自動的に更新されるように なっています． 3次元データ用の array_ptr_3d() もあります．

オブジェクト内部の配列データは常に1次元バッファで管理されていますから， 内部バッファのアドレスを 単純なポインタ変数に取得して，要素にアクセスする事もできます．

    float *arr0_ptr = arr0.array_ptr();
    size_t i, j, ii;
    ii = 0;
    for ( i=0 ; i < arr0.row_length() ; i++ ) {       /* Y */
        for ( j=0 ; j < arr0.col_length() ; j++ ) {   /* X */
            arr0_ptr[ii] = 100 + 10*i + j;
            ii ++;
        }
    }
    arr0.dprint();

当然ですが，ポインタ変数で要素にアクセスする場合は バッファの境界チェックは一切行なわれません．

○実習16○

上記の3つのコードで，配列の要素にどれも同じ値が入る事を 確認してください．

一部分だけ切り出してコピー，型はdoubleへ変換

上で作った float 型の2次元配列 arr0 の一部を切り出し， double 型の配列オブジェクトへコピーしてみましょう．

    mdarray_double arr1(false);
    arr1 = arr0.sectionf("1:4, *");
    arr1.dprint();

「"1:4, *"」が IDL 的な表記で， 「横方向は 1〜4 までの要素，縦方向はすべての要素」を示しています． もちろん，要素の番号は 0-indexed です． なお，FITSファイルのヘッダなどに，IRAF的表記で 1-indexed で配列の範囲が 書かれている事もあります．その場合は，"[…,…]" のように 書けば 1-indexed として扱われます． 今の場合は，"[2:5, *]" と書く事もできるわけです．

○実習17○

実習16のコードのつづきに上記のコードを追加し， 次のような値になる事を確認してください．

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 に貼り付けてみましょう．

    arr0.pastef(arr1, "*, 2:3");
    arr0.dprint();

pastef() の2つ目の引数で，arr0 のどの部分に 貼り付けるのかを指定します． 結果は次のとおり．左下の 4 × 2 の部分が貼り付けた部分です．

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 }
}

では，この状態からさらに同じ部分を arr1 で 割ってみましょう．

    arr0.dividef(arr1, "*, 2:3");
    arr0.dprint();

○実習18○

実習17のコードのつづきに上記のコード(ペーストおよび割り算)を追加し， 次のような値になる事を確認してください．

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 }
}

アレイに対して数学関数を使う

sli/mdarray_math.h を include すると， mdarray クラスやその継承クラスの配列オブジェクトに対し， 数学関数を使う事ができます． math.h で定義されているほとんどの数学関数をサポートしています． 例えば，配列の全要素を 2乗する場合は，次のように書く事ができます．

    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 }
}

統計値を得る

sli/mdarray_statistics.h を include すると， mdarray クラスやその継承クラスの配列オブジェクトに対しする 統計用の関数 を使う事ができます． すべての関数のコードは， sli/mdarray_statistics.h に記述されているため， SLLIB がインストールされている環境ならコードの確認・流用が手軽に行なえます．内容を確認してみましょう．

$ less /usr/local/include/sli/mdarray_statistics.h

arr0 の代表的な統計量を計算してみましょう:

    /* 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
}

y 方向で「本物 median」を得る

x方向，y方向あるいは z方向に統計をとり， 元の配列から次元数を1次元減らし，統計値を要素にセットした配列 を取得する事ができます．

例えば，これまでのコードの続きで arr0 を y 方向に median をとる場合は次のように書きます．

    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 }
    }

○実習19○

数学関数，統計値について示したコードの動作確認を行なってください．

パフォーマンス・テスト

SLLIB のメンバ関数や関数には，特別な最適化がなされているものもあります． そのパフォーマンスを確認してみましょう．

SLLIBのトランスポーズを試す

第一部の 「ハードウェアを意識しよう: メモリアクセス」 では，巨大二次元配列の x と y との入れ換えを行なう場合， 単純な方法ではパフォーマンスが出ない事を説明しました． SLLIB はこのトランスポーズを行なうための最適化されたコードを持っていますので， 試してみましょう．

以下のコードは，2次元配列 arr0 から， トランスポーズされた2次元配列 arr1 を作るものですが， 次の2つの場合について実行速度を測るものです．

• 2次元配列用のポインタ配列を使って素直に x と y とを 1ピクセルずつ入れ替えてコピーした場合
• mdarray クラスの transpose_xy_copy() を使った場合

#include <sli/stdstreamio.h>
#include <sli/mdarray.h>
#include <stdlib.h>
#include <time.h>
using namespace sli;

#define USE_SIMPLE_METHOD 1

const size_t Width = 4096*1;
const size_t Height = 4096*1;

int main( int argc, char *argv[] )
{
    stdstreamio sio;
    mdarray_double arr0, arr1;
    double tm0, tm1;
    size_t i;

    arr0.resize_2d(Width,Height);
    arr1.resize_2d(Width,Height);

    double *const *arr0_ptr2d = arr0.array_ptr_2d(true);
#ifdef USE_SIMPLE_METHOD
    double *const *arr1_ptr2d = arr1.array_ptr_2d(true);
#endif

    /* set random number */
    srand48(12345);
    for ( i=0 ; i < arr0.row_length() ; i++ ) {         /* Y */
        size_t j;
        for ( j=0 ; j < arr0.col_length() ; j++ ) {     /* X */
            arr0_ptr2d[i][j] = drand48();
        }
    }

    /* transpose: arr0 => arr1 */
    if ( arr0.col_length() == arr1.row_length() &&
         arr0.row_length() == arr1.col_length() ) {

#ifdef USE_SIMPLE_METHOD

        const size_t len_sq = arr0.col_length();
        tm0 = clock() / (double)CLOCKS_PER_SEC;
        for ( i=0 ; i < len_sq ; i++ ) {
            size_t j;
            for ( j=0 ; j < len_sq ; j++ ) {
                arr1_ptr2d[i][j] = arr0_ptr2d[j][i];
            }
        }
        tm1 = clock() / (double)CLOCKS_PER_SEC;
        sio.printf("direct 2d ptr transpose: %g\n",tm1 - tm0);

#else

        tm0 = clock() / (double)CLOCKS_PER_SEC;
        arr0.transpose_xy_copy(&arr1);
        tm1 = clock() / (double)CLOCKS_PER_SEC;
        sio.printf("transpose_xy_copy(): %g\n",tm1 - tm0);

#endif

    }

    arr0.dprint();
    arr1.dprint();

    return 0;
}

○実習20○

上記のコードについて，#define USE_SIMPLE_METHOD 1 を有効にした場合，無効にした場合の動作速度を確認してください．

SLLIBの md_median() を試す

SLLIB の 多次元配列用の統計用関数 によって計算される median は，IRAF のような近似値ではなく， 本物の median です． これも最適化がなされているので，その速度性能を計ってみましょう．

以下のコードは，次の2つの場合について実行速度を測るものです．

• std::sort() を使って median を求めた場合．
• mdarray_statistics.h の md_median() を使った場合

#include <sli/stdstreamio.h>
#include <sli/mdarray.h>
#include <sli/mdarray_statistics.h>
#include <stdlib.h>
#include <time.h>
#include <algorithm>                        /* std::sort() を使うために必要 */
using namespace sli;

const size_t Width = 4096*1;
const size_t Height = 4096*1;

int main( int argc, char *argv[] )
{
    stdstreamio sio;
    mdarray_double arr0;
    double tm0, tm1;
    size_t i;

    arr0.resize_2d(Width,Height);

    double *arr0_ptr = arr0.array_ptr();

    /* set random number */
    srand48(12345);
    for ( i=0 ; i < arr0.length() ; i++ ) {
        arr0_ptr[i] = drand48();
    }

    /* get median */
#define USE_STD_SORT 1

#ifdef USE_STD_SORT

    tm0 = clock() / (double)CLOCKS_PER_SEC;
    std::sort(arr0_ptr, arr0_ptr + arr0.length());
    double median = ( arr0[arr0.length() / 2 - 1] +
                      arr0[arr0.length() / 2] ) / 2.0;
    tm1 = clock() / (double)CLOCKS_PER_SEC;
    sio.printf("using std::sort(): %g\n",tm1 - tm0);
    sio.printf("median = %.15g\n",median);

#else

    tm0 = clock() / (double)CLOCKS_PER_SEC;
    double median = md_median(arr0);
    tm1 = clock() / (double)CLOCKS_PER_SEC;
    sio.printf("using md_median(): %g\n",tm1 - tm0);
    sio.printf("median = %.15g\n",median);

#endif

    return 0;
}

○実習21○

上記のコードについて，#define USE_STD_SORT 1 を有効にした場合，無効にした場合の動作速度を確認してください．

SFITSIO

いよいよ FITS ファイルを扱います． FITS ファイルの規約を理解していれば， LIBC あるいは SLLIB の関数を使って自力で I/O を行なう事もできますが， 最近は FITS の規約も複雑化している事もあり， FITS I/O のためのライブラリを使うのが一般的です．

FITS I/O のライブラリは，Python では PyFITS， IDL では AstroLib (MRDFITS)，C では CFITSIO，C++ では CCFITS などがあります． ここで扱う SFITSIO は PyFITS や IDL の MRDFITS に近く， オブジェクト指向の「オンメモリ型」の FITS I/O ライブラリです． ディスクベースの FITS I/O のための API も提供しており， ヘッダのみの高速読み取りのような応用が可能です． SLLIB と同様，プログラマが科学の本質部分のコーディングに集中できる 開発環境を提供する事を目的に開発され， JAXA 宇宙科学研究所の基幹部分でも使われています．

右図は SFITSIO を使った場合の典型的なプログラムが， どのような構成になるかを示しています． この図に示すとおり，SFITSIO は SLLIB の上に載っかっているライブラリです． この2つのライブラリは同じような設計思想を持ちますので， SLLIB を使った経験があれば SFITSIO の使用も簡単です．

IRAF の imcopy を作る

IRAF の imcopy タスクは，FITS ファイル(画像)の一部を読み取り， 別の FITS ファイルへ保存するものでした．

これと同じようなものを作ってみましょう． なんだか難しいような気がしますが， SFITSIO では「読んで」「書く」， 2つのメンバ関数 を使うだけです．

#include <sli/stdstreamio.h>
#include <sli/fitscc.h>       /* SFITSIO 使用時はこれを includeすればOK */
using namespace sli;

int main( int argc, char *argv[] )
{
    int return_status = -1;
    stdstreamio sio;                                    /* standard I/O */
    fitscc fits;                                        /* FITS object  */

    /* FITSファイルを読み，fits オブジェクトに内容を保持させる */
    if ( 1 < argc ) {
        const char *in_file = argv[1];
        if ( fits.read_stream(in_file) < 0 ) {
            sio.eprintf("[ERROR] fits.read_stream(\"%s\") failed\n",in_file);
            goto quit;
        }
    }

    /* fits オブジェクトの内容を FITSファイルとして保存する */
    if ( 2 < argc ) {
        const char *out_file = argv[2];
        if ( fits.write_stream(out_file) < 0 ) {
            sio.eprintf("[ERROR] fits.write_stream(\"%s\") failed\n",out_file);
            goto quit;
        }
    }

    return_status = 0;
 quit:
    return return_status;
}

上のコードはエラーのチェックが入っているのでやや長くなっていますが， 本質的な部分は

    fits.read_stream(argv[1]);
    fits.write_stream(argv[2]);

の2行だけです．

SFITSIO を使ったコードをコンパイルする時は，次のように 「-lsfitsio」を最後の引数として追加します．

$ s++ read_and_write.cc -lsfitsio

これで imcopy に相当するものができました．

○実習22○

1. 上で示したコードをコンパイルし， 次のように画像データ(Image HDU)を含む FITS ファイル M83_DSS2_R_A2EE.fits.gz を読み込み， 別のファイルへ保存し，ds9，fv で FITS の内容を 確認してください．
   • すべてをコピー

     $ ./read_and_write M83_DSS2_R_A2EE.fits.gz out_0.fits
     

   • Primary HDU のみコピー

     $ ./read_and_write M83_DSS2_R_A2EE.fits.gz'[0]' out_1.fits
     

     この場合，続く HDU へはアクセスせず， Primary HDU の読み取りが完了した時点でストリームをクローズします． したがって，続く HDU が破損していた場合や ゴミデータが続くような不正な FITS ファイル(一時期のCOMICSのFITSファイルが有名)でも，[0] を指定していれば問題なく読む事ができます．
   • Primary HDU の画像の一部領域をコピー

     $ ./read_and_write M83_DSS2_R_A2EE.fits.gz'[1:100,*]' out_2.fits
     

     厳密には次のように書く:

     $ ./read_and_write M83_DSS2_R_A2EE.fits.gz'[0[1:100,*]]' out_3.fits
     

     0-indexed で指定する場合は丸括弧を使う:

     $ ./read_and_write M83_DSS2_R_A2EE.fits.gz'[0(0:99,*)]' out_4.fits
     

     画像データの領域選択については，次元数の制限はありません． 例えば，[0[*,*,1]] と指定すれば， 3次元データの最初のレイヤだけを選択する事ができます．
2. 次のようにバイナリテーブルを含む FITS ファイル tsd4_sample.fits.gz を読み込み， 別のファイルへ保存し，fv で FITS の内容を確認してください． 下の図のように，fv の GUI で各 HDU の「All」をクリックすると， テーブルの内容を確認する事ができます．
   
   • すべてをコピー

     $ ./read_and_write tsd4_sample.fits.gz out_a.fits
     

   • Primary HDUと "FIS_HK" HDU のみ取り出す．

     $ ./read_and_write tsd4_sample.fits.gz'[0; FIS_HK]' out_b.fits
     

   • 最初の3つのHDUだけを取り出す

     $ ./read_and_write tsd4_sample.fits.gz'[*{3}]' out_c.fits
     

   • HDU名に文字列 "HK" を含む HDU だけを取り出す．

     $ ./read_and_write tsd4_sample.fits.gz'[*HK*]' out_d.fits
     

   • すべての HDU について，テーブルのカラム名に文字列 "TI" を含むカラムと カラム "DET" のみを取り出す．

     $ ./read_and_write tsd4_sample.fits.gz'[*[*TI*;DET]]' out_e.fits
     

   • 上記の条件からさらに，最初の 100 行のみ取り出す．

     $ ./read_and_write tsd4_sample.fits.gz'[*[*TI*;DET, 1:100]]' out_f.fits
     

     0-indexed で指定する場合は次のとおり．

     $ ./read_and_write tsd4_sample.fits.gz'[*(*TI*;DET, 0:99)]' out_g.fits
     

参考:
ファイル名の後の […] を使った FITS ファイルの部分読み出し は，FITS ファイルの内容すべてをメモリに取り込む事はなく， ディスクのシーケンシャルリードと順方向のシーク(可能な場合) により実現しています． したがって， メモリに載らないような巨大な FITS ファイルでも(圧縮ファイルであっても)， 必要な部分を読み出す事ができます．

fitscc オブジェクトの内容と目的のデータへのアクセス

    fitscc fits;
    fits.read_stream(in_file);

さきほど出てきたこのコードは，FITS ファイルの内容のすべて ([…] をつけた場合は内容の一部) を読み取り，fits というオブジェクト(が管理するメモリ) に取り込むものです．

SLLIB の文字列配列 (tarray_tstringクラス)の扱いでは， 文字列配列オブジェクトの内部で複数の tstring オブジェクト が管理されており，それを引き出すのに array[i] のように した事を思い出してください． SFITSIO の場合もこれと同様な作りになっていて， fitsccクラス のオブジェクトの場合は， 内部で管理しているオブジェクトが入れ子になっていて， 複数のメンバ関数を使って必要な内部オブジェクトを引き出します．

次の図は，fitscc オブジェクト内部についての模式図で， FITS ファイルの内容のヘッダユニット・データユニットの内容が， 対応するクラスのオブジェクトで管理されています． 要するに， FITSファイルの内容がオブジェクトの構造に転写されている わけです．

次の表に，SFITSIO で使うクラスについての一覧を示します．

FITS のデータ要素へのアクセスの例を示します:

    printf("TELESCOP = %s\n",
           fits.image("Primary").header("TELESCOP").svalue());

この場合，image("Primary") は fits_image クラス のオブジェクトを引き出し，header("TELESCOP") は fits_header_record クラス のオブジェクトを引き出しています． 最後の svalue() で ヘッダ TELESCOP の文字列の値を取り出しています．

多くのヘッダレコードにアクセスする場合などには， fits.image(... と書くと長いしパフォーマンスにも 良くありません．この場合は，参照を使って短く書く事ができます．

    fits_image &img = fits.image("Primary");

    printf("TELESCOP = %s\n", img.header("TELESCOP").svalue());
    printf("GAIN = %g\n", img.header("GAIN").dvalue());

FITS ヘッダの値や，画像データのピクセル値などの， プリミティブなデータへのアクセスも， パターンが決まっているので，表にまとめておきます．

以上で FITS のほとんどのデータ要素へアクセスできるようになります． 代表的なものをまとめておきましょう．

なお，dvalue()，assign() 等を使う方法では， Image HDU の BZERO値，BSCALE値，BLANK値， Binary Table HDU または ASCII Table HDU の ZEROn値，TSCALn値，TNULLn値 が考慮され，ライブラリ側で変換処理が行なわれます (lvalue()の場合は四捨五入された値を得ます)．

SFITSIO で定義される定数と型

SFITSIO では， 「fitscc オブジェクトの内容と目的のデータへのアクセス」 でみてきたクラスのほか， 型や型を示す整数値などの定数が /usr/local/include/sli/fits.h で定義されています． ここでは，それらについてまとめておきます．

なお，これらの定数や型は，すべて namespace sli の中で定義されています． 例えば，「FITS::IMAGE_HDU」は完全な表記では 「sli::FITS::IMAGE_HDU」です．

• HDUのタイプを示す定数

  型	定数	実際の値	意味
  int	FITS::ANY_HDU	127	不明なタイプのHDU
  int	FITS::IMAGE_HDU	0	Image HDU
  int	FITS::BINARY_TABLE_HDU	2	Binary Table HDU
  int	FITS::ASCII_TABLE_HDU	1	ASCII Table HDU

• 長さを示す定数

  型	定数	実際の値	意味
  long	FITS::ALL	longの最大値	すべての配列要素
  long	FITS::INDEF	longの最小値	未定義を示す

• ヘッダレコードの状態を示す定数

  型	定数	実際の値	意味
  int	FITS::NULL_RECORD	0	すべてが空白文字である
  int	FITS::NORMAL_RECORD	1	「キーワード = 値 / コメント」の形式を示す
  int	FITS::DESCRIPTION_RECORD	2	「キーワード コメント」の形式を示す

• データ型を示す定数

  型	定数	実際の値	意味	ヘッダ	Image	Binary Table
  int	FITS::ANY_T	127	不明な型
  int	FITS::ASCII_T	65('A')	FITSの文字型(char型)	Yes	-	Yes
  int	FITS::BIT_T	88('X')	FITSのbit型	-	-	Yes
  int	FITS::LOGICAL_T	76('L')	FITSの論理値型(unsigned char型)	Yes	-	Yes
  int	FITS::BYTE_T	66('B')	FITSの符号なし1バイト整数型(unsigned char型)	-	Yes	Yes
  int	FITS::SHORT_T	73('I')	FITSの2バイト整数型(16-bit整数型)	-	Yes	Yes
  int	FITS::LONG_T	74('J')	FITSの4バイト整数型(32-bit整数型)	-	Yes	Yes
  int	FITS::LONGLONG_T	75('K')	FITSの8バイト整数型(64-bit整数型)	Yes	Yes	Yes
  int	FITS::FLOAT_T	69('E')	float型	-	Yes	Yes
  int	FITS::DOUBLE_T	68('D')	double型	Yes	Yes	Yes
  int	FITS::COMPLEX_T	67('C')	FITSの複素数型(float _Complex型)	-	-	Yes
  int	FITS::DOUBLECOMPLEX_T	77('M')	FITSの倍精度複素数型(double _Complex型)	Yes	-	Yes
  int	FITS::LONGARRDESC_T	80('P')	32-bit配列記述子(32-bit整数×2)	-	-	Yes
  int	FITS::LLONGARRDESC_T	81('Q')	64-bit配列記述子(64-bit整数×2)	-	-	Yes

• SFITSIOで定義されるデータ型

  SFITSIOで使う型	実際の型	Image	Binary Table
  fits::ascii_t	char	-	Yes
  fits::bit_t	struct {   unsigned char elem7:1;   unsigned char elem6:1;   unsigned char elem5:1;   unsigned char elem4:1;   unsigned char elem3:1;   unsigned char elem2:1;   unsigned char elem1:1;   unsigned char elem0:1; };	-	Yes
  fits::logical_t	uint8_t	-	Yes
  fits::byte_t	uint8_t	Yes	Yes
  fits::short_t	int16_t	Yes	Yes
  fits::long_t	int32_t	Yes	Yes
  fits::longlong_t	int64_t	Yes	Yes
  fits::float_t	float	Yes	Yes
  fits::double_t	double	Yes	Yes
  fits::complex_t	float _Complex	-	Yes
  fits::doublecomplex_t	double _Complex	-	Yes
  fits::longarrdesc_t	struct {   uint32_t length;   uint32_t offset; };	-	Yes
  fits::llongarrdesc_t	struct {   uint64_t length;   uint64_t offset; };	-	Yes

• ヘッダレコードを定義するための型(構造体)

  SFITSIOで使う型	実際の型
  fits::header_def	struct {   const char *keyword;   const char *value;   const char *comment; };

• バイナリテーブル・ASCIIテーブルを定義するための型(構造体)

  SFITSIOで使う型	実際の型
  fits::table_def	struct {   const char *ttype;   const char *ttype_comment;   const char *talas;   const char *telem;   const char *tunit;   const char *tunit_comment;   const char *tdisp;   const char *tform;   const char *tdim;   const char *tnull;   const char *tzero;   const char *tscal; };

新規 FITS ファイルの作成

では，実際に新しい FITS ファイルを作ってみましょう．

#include <stdio.h>
#include <sli/fitscc.h>        /* required by every source that uses SFITSIO */
using namespace sli;

int main( int argc, char *argv[] )
{
    fitscc fits;                 /* fitscc object that expresses a FITS file */
    const long width = 300, height = 300;                   /* size of image */
    long len_written, i;

    /* fits オブジェクトに Primary HDU を追加 */
    fits.append_image("Primary",0, FITS::SHORT_T, width,height);
    fits_image &pri = fits.image("Primary");

    /* ヘッダに "EXPOSURE" キーワードを追加 */
    pri.header("EXPOSURE").assign(1500.).assign_comment("Total Exposure Time");

    /* BZERO，BSCALE の設定と初期化 */
    pri.assign_bzero(32768.0).assign_bscale(1.0).fill(0);
    pri.assign_blank(65535);
    pri.assign_bunit("ADU");

    /* ピクセル値をセット */
    for ( i=0 ; i < height ; i++ ) {
        long j;
        for ( j=0 ; j < width ; j++ ) pri.assign(j + i, j,i);
    }

    /* FITS ファイルを出力 */
    len_written = fits.write_stream("testfile.fits");
    printf("saved %ld bytes.\n", len_written);

    if ( len_written < 0 ) return -1;
    else return 0;
}

新規の Image HDU を作るには， fitsccクラス のオブジェクトを作って，

    fits.append_image(HDU名,HDUバージョン, データ型, xの長さ,yの長さ,zの長さ);

とするだけで，あとは fitscc オブジェクトの内容と目的のデータへのアクセス で示したようにメンバ関数を使ってアクセスするだけです． ヘッダへのアクセスや，画像のピクセル値への高速アクセス などについては，この後で紹介します．

○実習23○

1. 上のコードをコンパイルし，新規 FITS を作ってみてください．
2. BZERO，BSCALE の設定を削除し， /usr/local/include/sli/fits.h をみて，4バイト実数(float)型の新規データを作成するよう コードを改造してください．
3. 「fits.append_image(…);」の直後に， 次のコードを挟んでコンパイルし，動作確認を行なってください．

       fits::header_def defs[] = { {"TELESCOP", "'CASSIOPEIA'", "Telescope name"},
                                   {"OBSERVAT", "'NAOJ'", "Observatory name"},
                                   {"RA",       "", "[deg] Target position"},
                                   {"DEC",      "", "[deg] Target position"},
                                   {"COMMENT",  "-------------------------------"},
                                   {NULL} };
       /* ヘッダの初期化 */
       fits.image("Primary").header_init(defs);
   

画像データへのアクセス

画素の大きさは， fits_imageクラスの col_length()， row_length()， layer_length() で調べる事ができます． 例えば，画素の大きさを表示するには，次のように書けます．

    fits_image &primary = fits.image("Primary");
    printf("カラムの数 : %ld\n",primary.col_length());
    printf("ロウの数   : %ld\n",primary.row_length());
    printf("レイヤの数 : %ld\n",primary.layer_length());

読み出しには， fits_imageクラスの dvalue()， lvalue()， llvalue() のいずれかを使い， 書き込みには assign() を使います． これらは，ヘッダの BZERO，BSCALE の値による変換もやってくれます．

dvalue() を使うと， イメージデータの型にかかわらず，doubleで読む事ができます． 座標値は 0-indexed です．

    double pixel_val;
    pixel_val = primary.dvalue(x,y,z);   /* 読み */
    primary.assign(pixel_val,x1,y1,z1);  /* 書き */

もちろん，整数値で読み書きする事もできます． lvalue()，llvalue() でそれぞれ long，long long を返してくれます．

    long pixel_val;
    pixel_val = primary.lvalue(x,y,z);   /* 読み */
    primary.assign(pixel_val,x1,y1,z1);  /* 書き */

dvalue() や assign() は高いレベルのAPIですので安全性は 高いですが，その分呼び出しのオーバヘッドは小さくありません． より高速なアクセスについてはこの後で紹介します．

画像データの型変換と高速アクセス

オブジェクト内部データへのポインタ変数を使った高速アクセスも可能です． ただし， fits_imageクラス の resize_2d()メンバ関数 を使った場合などで画素のサイズが変更になった場合は，データ配列のアドレスが変わ るので注意が必要です．

高速アクセスのためには，値を読む時の BZERO値とBSCALE値による変換がいらないように 内部データを変換してしまうと良いでしょう． 変換には， fits_imageクラス の convert_type()メンバ関数 を使います．例えば，次のようにします．

    fits_image &primary = fits.image("Primary");
    primary.convert_type(FITS::FLOAT_T);

以上で，どの形式のデータもBZERO値が0，BSCALE値が1のfloat型に変換されます． この後， fits_imageクラス の float_t_ptr_2d()メンバ関数，float_t_ptr()メンバ関数等 を使って，データのアドレスを得ます(それぞれの型に対応したメンバ関数があります)．

    fits::float_t *const *ptr_2d;
    ptr_2d = primary.float_t_ptr_2d(true);

上記の場合は「ptr_2d[y][x]」の形で高速アクセスが可能です．

    fits::float_t *ptr;
    ptr = primary.float_t_ptr();

この場合は 「ptr[x + y * col_length]」の形でアクセスする事になります．

◎課題4◎

1. 任意の FITS ファイルを読み込み，Primary HDU の画像の型を 4バイト実数に変換し， ピクセル値の平均値を求めるプログラムを作ってください． a0117.fits.gz を使って動作確認を行なってください．

ヘッダへのアクセスの基本

ヘッダへのアクセスへの基本はこれまでみてきたとおりです:

    printf("TELESCOP = %s\n", fits.hdu("Primary").header("TELESCOP").svalue());

と書きます．HDUの指定は上記のように名前で指定するか， 数値(プライマリHDUの場合は「.hdu(0L)」)でもかまいません． なお，この例では「fits.hdu(…)」と書いていますが， これは HDU のタイプが不明な時の書き方です． HDU のタイプがわかっている時は， 「fits.image(…)」と書いてもかまいません．

値の読み出しには dvalue()，lvalue()，llvalue()， bvalue()，svalue() のいずれかを使います． それぞれ，double，long， long long，bool，const char * の各型に対応します．

逆に，書き込む場合は，新規ヘッダレコード追加の場合も 値やコメントの更新の場合も， assign() や assign_comment() を使います．

    fits.hdu("Primary").header("TELESCOP").assign("HST")
                                          .assign_comment("Telescope Name");
    fits.hdu("Primary").headerf("OBJECT%d",n).assignf("%s-%d",obj[i],j)
                       .assignf_comment("Name of the object No.%d",n);

のように書きます．LIBC の printf() 関数の書き方が いたるところで使えるので，大変便利です． HISTORY などの記述レコードの追加も簡単で， header_append() メンバ関数に2つの引数を与え， このように書きます:

    fits.hdu("Primary").header_append("HISTORY","step-0: done.");

ヘッダのキーワードが存在し，値が存在するかどうかをチェックするには， header_value_length()メンバ関数が便利です．

    if ( 0 < fits.hdu("Primary").header_value_length("TELESCOP") ) {
        /* OK */
    }

このメンバ関数は，記述レコード以外から与えられたキーワードを検索し， 存在しない場合は負値，存在する場合はそのヘッダレコードの値の長さ (文字列では「'」も含む)を返します．

ヘッダのキーワードと値が存在する事がわかったら， 値がどの型にあてはまるかをチェックしたい事もあります． その場合は，type() を使います．

    int r_type = fits.hdu("Primary").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::LOGICAL_T ) {
        /* 論理値 */
    } else if ( r_type == FITS::ASCII_T ) {
        /* 文字列 */
    } else {
        /* 予期せぬエラー */
    }

なお，Primary HDUの場合は必ずイメージを扱うHDUですので， fits.hdu(...) は fits.image(...) としてもかまいません．

s++ コマンドで表示されるのサンプルコードもご覧ください．

ヘッダの編集

ヘッダの初期化(header_init())， ヘッダレコードの追加(header_append_records())， 挿入(header_insert_records())， 削除(header_erase_records()) を行なうメンバ関数が fits_hduクラス で用意されています (fits_imageクラス や fits_tableクラス は fits_hduクラス を継承しているので，同様に使えます)．

次のコードは，一方のFITSのヘッダレコードのすべての内容を， もう1方にコピーする例です．

    fits_out.image("Primary")
            .header_append_records( fits_in.image("Primary").header() );

このように，引数が無い形 .header() は，ヘッダ全体を表現しており，同様にして， header_init()メンバ関数， header_insert_records()メンバ関数 に与える事ができます．

ヘッダレコードのすべての内容ではなく， 1つのヘッダレコードをコピーする場合は，次のようにします．

    fits_out.image("Primary")
            .header_append( fits_in.image("Primary").header("TELESCOP") );

以下に，既存の FITS ファイルのヘッダを，新規 FITS ファイルに コピーするコードを示します． 「IRAF の imcopy を作る」 のコードを少し手直ししただけです．

#include <sli/stdstreamio.h>
#include <sli/fitscc.h>       /* SFITSIO 使用時はこれを includeすればOK */
using namespace sli;

int main( int argc, char *argv[] )
{
    int return_status = -1;
    stdstreamio sio;                                /* standard I/O */
    fitscc in_fits;                                 /* FITS object (in) */

    /* FITS ファイルを読む */
    if ( 1 < argc ) {
        const char *in_file = argv[1];
        if ( in_fits.read_stream(in_file) < 0 ) {
            sio.eprintf("[ERROR] in_fits.read_stream(\"%s\") failed\n",in_file);
            goto quit;
        }
    }

    /* ヘッダをコピーした新規 FITSファイルを作成 */
    if ( 2 < argc ) {
        const char *out_file = argv[2];
        fitscc out_fits;                            /* FITS object (out) */

        const fits_image &in_img = in_fits.image("Primary");

        out_fits.append_image(in_img.extname(), in_img.extver(), FITS::FLOAT_T,
                              in_img.col_length(), in_img.row_length());
        out_fits.image("Primary")
                .header_append_records( in_img.header(), false );

        if ( out_fits.write_stream(out_file) < 0 ) {
            sio.eprintf("[ERROR] out_fits.write_stream(\"%s\") failed\n",out_file);
            goto quit;
        }
    }

    return_status = 0;
 quit:
    return return_status;
}

オブジェクト指向の特徴は，このように それぞれのオブジェクトが FITS の部品を表現していて， メンバ関数を使って部品単位で編集ができるという事です．

○実習24○

上記のコードをコンパイルし，動作を確認してください． 元の FITS ファイルは M83_DSS2_R_A2EE.fits.gz を利用してください． なお，自身のものと第1引数で指定したものとで ヘッダキーワードが重複するレコードがある場合，第1引数のものは 無視されます．header_append_records() の2番目の引数は， この重複が起きた場合に警告するかどうかを指定しています (falseだと報告しない)．

画像データのコピー & ペースト

fits_imageクラス には IDL/Python風の表記を使って矩形の領域のコピーとペーストを行なうためのメンバ関数， copyf() と pastef() が用意されています． 次に示すのは，座標(0,0)から(99,99)までの矩形領域を， 座標(100,100)へコピーする例です．

    fits_image &primary = fits.image("Primary");
    fits_image copy_buf;
    primary.copyf(&copy_buf, "0:99, 0:99");
    primary.pastef(copy_buf, "100:*, 100:*");

プログラマが作った fits_image のオブジェクト copy_buf を コピーバッファとして使っています． その copy_buf オブジェクトをメンバ関数に与えて， コピー & ペーストを行ないます． コピーバッファはいくつでも持つ事ができますし， オブジェクト間のコピー & ペーストも簡単です．

pastef() のかわりに addf()，subtractf()， multiplyf()，dividef() を使うと， 足し算，引き算，掛け算，割り算もできます．

本格的な画像解析ツールを開発するためのヒント

速度が重視される本格的な画像解析ツールを開発する場合には， SLLIB の単純な実装によるAPIを使ってコーディングを行なう事で， 全体的なパフォーマンスアップが可能です．

次に SLLIB ベースの解析を行なう手順を述べます． 考え方としては，IDL や Python + numpy の場合と同様， FITSの画像データを単純な配列オブジェクトとして プロセッシングを行なうというものです．

まず，FITSファイルを読み，データ型をfloat型またはdouble型に変換します．

    fitscc fits;
    fits.read_stream("largefile.fits");
    fits_image &prim = fits.image(0L);
    prim.convert_type(FITS::FLOAT_T);

次に， mdarray_floatクラス で参照を取得します．

    mdarray_float &prim_arr = prim.float_array();

この後，配列に対しては SLLIB の mdarray クラスの API を使い， FITSヘッダに対しては SFITSIO の API を使います． 例えば，ゲインをかけて画像の一部の 統計をとるには 次のように書けます．

    double gain = prim.header("GAIN").dvalue();
    prim_arr *= gain;

    double meanabsdev, stddev;
    /* get mean, variance, skewness, kurtosis, and stddev */
    mdarray_double moment = md_moment(prim_arr.sectionf("2:21, *"), false,
                                      &meanabsdev, &stddev);
    printf("mean = %g\n", moment[0]);
    printf("variance = %g\n", moment[1]);
    printf("skewness = %g\n", moment[2]);
    printf("kurtosis = %g\n", moment[3]);
    printf("meanabsdev = %g\n", meanabsdev);
    printf("stddev = %g\n", stddev);

なお，SLLIB は FITS とは無関係なライブラリなので， プリミティブ型の名称が SFITSIO の場合とは異なります． 次に対応関係を示します．

◎課題5◎

1. FITS ファイル a0117.fits.gz を読み込み，Primary HDU の画像について，全体，オーバスキャン領域の統計値(平均値，標準偏差)を md_moment()関数 を使って求めるプログラムを作ってください． オーバスキャン領域はヘッダの BIASSEC の値に記録されています． なお，統計値取得前にゲインをかける必要はありません．

WCSToolsのlibwcsとの連携

WCSTools (http://tdc-www.harvard.edu/wcstools/) は WCS (World Coordinate System) をはじめとした， 天文学のデータを効率良く扱うために開発されたC言語ライブラリです． FITSファイルや天体カタログファイルにも対応しており， SAO の Douglas J. Mink 氏が開発・メンテナンスしています．

ここでは，SFITSIOをlibwcsと組み合わせると， 簡単にWCSを扱える事を示してみます．

コードの例を示します．始めに，引数で指定された FITSファイルを読み込み， wcsinitn()関数を使って， プライマリHDUについて，FITSヘッダから WorldCoor構造体の オブジェクトを wcsに作成しています． この時に，プライマリHDUのすべてのヘッダレコードを1つの文字列として返して header_formatted_string()メンバ関数 を利用しています． 次の iswcs() 関数では， WorldCoor構造体のオブジェクトが正しくセットされたか どうかをチェックしています (iswcs()関数は，正しくセットされた場合は1，そうでない場合は 0を返します)． その後，pix2wcs()関数を使ってピクセル座標(0,0)について 世界座標(lon, lat)を求め， さらに wcs2pix() 関数で逆の変換をしています． 最後に当該ピクセルの値を表示し，wcs が使っているメモリを開放しています．

#include <sli/stdstreamio.h>
#include <sli/fitscc.h>
#include <math.h>
#include <libwcs/wcs.h>
using namespace sli;

int main( int argc, char *argv[] )
{
    int ret_status = -1;
    stdstreamio sio;
    fitscc fits;                                        /* FITS object */
    struct WorldCoor *wcs = NULL;                       /* WCS structure */
    double lon,lat, x,y, v;
    int off;

    if ( 1 < argc ) {

        const char *in_file = argv[1];

        /* FITSファイルを読む */
        if ( fits.read_stream(in_file) < 0 ) {
            sio.eprintf("[ERROR] fits.read_stream(\"%s\") failed\n",in_file);
            goto quit;
        }

        /* Primary HDU を示すオブジェクトへの参照の作成 */
        fits_image &pri = fits.image("Primary");

        /* wcs 構造体のオブジェクトを取得する */
        wcs = wcsinitn(pri.header_formatted_string(), NULL);

        if ( ! iswcs(wcs) ) {
            sio.eprintf("[ERROR] wcsinitn() failed\n");
            goto quit;
        }

        /* ピクセル座標 -> wcs の変換 */
        x = 1.0;  y = 1.0;
        pix2wcs(wcs, x, y, &lon, &lat);
        sio.printf("ra=%.8f dec=%.8f\n",lon,lat);

        /* wcs -> ピクセル座標 の変換 */
        wcs2pix(wcs, lon, lat, &x, &y, &off);
        sio.printf("x=%.8f y=%.8f\n",x,y);

        /* ピクセル値を読む (1-indexed) */
        v = pri.dvalue((long)floor(x-0.5),(long)floor(y-0.5));
        sio.printf("value=%.15g\n",v);

    }

    ret_status = 0;
quit:
    /* wcs 構造体のオブジェクトを開放 */
    if ( wcs != NULL ) wcsfree(wcs);
    return ret_status;
}

WCSToolsのこれ以上の情報については，ヘッダファイル wcs.h 等を参照してください． WCSToolsのヘッダファイルには，かなりしっかりとコメントが記載されており， それらを読むだけでもかなりの事ができるようになるはずです．

○実習25○

次の手順で，上記のコードを動かしてみてください．

1. 次の手順で wcstools-3.9.2.tar.gz を展開し，libwcs をビルドしてください．

   $ tar zxvf wcstools-3.9.2.tar.gz
   $ cd wcstools-3.9.2/libwcs
   $ make CFLAGS="-O -Wall"
   $ cd ../..
   

2. 上記のコードを次のコマンドでコンパイルし， 動作確認を行なってください．FITS ファイルは M83_DSS2_R_A2EE.fits.gz を使ってください．

   $ s++ -Iwcstools-3.9.2 wcs_sample.cc wcstools-3.9.2/libwcs/libwcs.a -lsfitsio
   

3. WCSを計算する座標を変更し，ds9 と同じ座標が得られている事を確認してください．

バイナリテーブルの作成とアクセス

「fitscc オブジェクトの内容と目的のデータへのアクセス」 で示したように，バイナリテーブルへのアクセスも 「なんとか.なんとか.なんとか」の形をとります:

    double value = fits.table("MY_TABLE").col("RA").dvalue(row_index);
    fits.table("MY_TABLE").col("RA").assign(value, row_index);

新規のバイナリテーブルを作成する場合は，次のように 構造体を使う方法があります:

    fitscc fits;
    const fits::table_def def[] = {
      /* ttype,comment,         talas,telem,tunit,comment,                */
      /*                                            tdisp,  tform, tdim, tnull */
      { "TIME","satellite time",   "",  "", "s","", "F16.3", "1D", "" },
      { "STATUS","status",         "",  "", "","",  "",      "8J", "", "-99999" },
      { "NAME","",                 "",  "", "","",  "",  "128A16", "(4,2)" },
      { NULL }
    };
    /* バイナリテーブルの作成(HDU名は"EVENT") */
    fits.append_table("EVENT",0, def);

しかし，SFITSIO は「FITSテンプレート」という機能を持ち， それを使う方が様々な応用が効くのでお勧めです． FITSテンプレートとは， FITSヘッダに類似した曖昧な文法でFITSの内容を定義できるテキスト形式のファイルで， これをもとにデータが入っていない新規 FITSファイル(オブジェクト)を 作成する事ができます． 例えば，上記の構造体で定義されたバイナリテーブルを， FITSテンプレートで書くと次のようになります．

XTENSION = BINTABLE
EXTNAME = EVENT
 TTYPE#  = 'TIME' / satellite time
 TFORM#  = '1D'
 TUNIT#  = 's'
 TDISP#  = 'F16.3'
 TTYPE#  = STATUS / status
 TFORM#  = '8J'
 TNULL#  = -99999
 TTYPE#  = 'NAME'
 TFORM#  = '128A16'
 TDIM#   = '(4,2)'

一方，コードは read_stream() のかわりに， read_template() を使うだけです．

#include <sli/stdstreamio.h>
#include <sli/fitscc.h>
using namespace sli;

int main( int argc, char *argv[] )
{ 
    int return_status = -1;
    stdstreamio sio;
    fitscc fits;

    if ( 1 < argc ) {
        const char *templ_filename = argv[1];

        /* FITS テンプレートを読み取る */
        if ( fits.read_template(templ_filename) < 0 ) {
            sio.eprintf("[ERROR] fits.read_template() failed\n");
            goto quit;
        }

        if ( fits.hdutype(1L) == FITS::BINARY_TABLE_HDU ) {
            long i, j;

            /* テーブルオブジェクトへの参照を取得 */
            fits_table &tbl = fits.table(1L);

            /* カラムオブジェクトへの参照を取得 */
            fits_table_col &col_time = tbl.col("TIME");
            fits_table_col &col_name = tbl.col("NAME");

            /* すべてのセルを NULL 値に初期化 */
            for ( i=0 ; i < tbl.col_length() ; i++ ) {
                tbl.col(i).assign_default(NAN);
            }

            /* 行を確保 */
            tbl.resize_rows(64);

#if 0
            /* セルに値をセット */
            for ( i=0 ; i < tbl.row_length() ; i++ ) {
                col_time.assign(1000 + i, i);
                for ( j=0 ; j < col_name.elem_length() ; j++ ) {
                    col_name.assign("None", i, j);
                }
            }
#endif
        }

    }

    if ( 2 < argc ) {
        const char *out_filename = argv[2];
        sio.printf("# writing ...\n");
        if ( fits.write_stream(out_filename) < 0 ) {
            sio.eprintf("[ERROR] fits.write_stream() failed\n");
            goto quit;
        }
    }

    return_status = 0;
 quit:
    return return_status;
}

○実習26○

1. 上記のFITSテンプレートとコードとをファイルに保存し， 上記コードを動作させて FITSテンプレートからFITSファイルを作ってみてください． 動作させる時は，1つめの引数に読み込ませるFITSテンプレートのファイルを， 2つめの引数に作成する FITS ファイルの名前を指定します． fv で作られた FITS ファイルを開き，ヘッダにコメントが自動で入っている事， バイナリテーブルのすべてのセルが NULL である事を確認してください．
2. 「#if 0」を「#if 1」にして動作させ， 値がセットされている事を確認してください．

ディスクベースのFITS I/O

これまでみてきたように，SFITSIO はオンメモリ型の FITS I/O が基本です． また，ファイル名直後の […] の表現でディスクベースでの "部分読み"もできるので，ゴミがついているファイルや 巨大なファイルでも読み取りは可能です．

しかし，場合によっては，ヘッダのみを高速に読み取りたい事や， ディスク上の画像データにランダムにアクセスしたい事もあります． SFITSIO ではそのような事も想定し， SLLIB でオープンしたストリームに対して， ヘッダユニット，データユニットにアクセスするための API を 提供しています．

最も簡単な例として，ヘッダ部分のみ読むコードを示します． このコードは，ヘッダの部分だけをストリームから読み，ストリームを閉じるので， たとえ圧縮ファイルであっても， ネットワーク上のファイルであっても， 非常に高速に動作します． ポイントは fits_headerクラス のメンバ関数を利用する事です．

#include <sli/stdstreamio.h>
#include <sli/digeststreamio.h>
#include <sli/fitscc.h>
using namespace sli;

int main( int argc, char *argv[] )
{
    int ret_status = -1;
    stdstreamio sio;
    digeststreamio f_in;

    if ( 1 < argc ) {

        const char *in_file = argv[1];
        fits_header hdr;                        /* ヘッダオブジェクト */

        /* ストリームをオープン */
        if ( f_in.open("r", in_file) < 0 ) {
            sio.eprintf("[ERROR] cannot open file: %s\n",in_file);
            goto quit;
        }

        /* Primary HDU のヘッダを読む */
        hdr.read_stream(f_in);

        /* ヘッダの内容を表示 */
        sio.printf("BITPIX = %ld\n", hdr.at("BITPIX").lvalue());

        long naxis = hdr.at("NAXIS").lvalue();
        long i;
        for ( i=0 ; i < naxis ; i++ ) {
            tstring key;
            key.printf("NAXIS%ld",i+1);
            sio.printf("%s = %ld\n", key.cstr(), hdr.at(key.cstr()).lvalue());
        }

        sio.printf("TELESCOP = %s\n", hdr.at("TELESCOP").svalue());
        sio.printf("DATE-OBS = %s\n", hdr.at("DATE-OBS").svalue());
        sio.printf("\n");

        sio.printf("[all string]\n");
        sio.printf("%s\n", hdr.formatted_string());

    }

    ret_status = 0;
quit:
    return ret_status;
}

「hdr.read_stream(f_in);」 の時点でファイルポインタは次のData Unitの先頭にあるので， f_in.read(...) などで Data Unit にアクセスが可能です． その場合， fits_headerクラス のオブジェクト hdr は Primary HDUの全ヘッダの内容を保持していますから， NAXIS1，NAXIS2 等の値の取得も簡単です．

Data Unit の読み飛ばしも簡単です (ただし，ここは圧縮ファイルだと速度は出ません):

        hdr.skip_data_stream(f_in);

このようにすると，hdr オブジェクトに入っている ヘッダ情報に基づいて Data Unit 部分をスキップします．

次の HDU のヘッダも同様に読み取りができます．

        hdr.read_stream(f_in);

必要なものを読んだら，いつでもストリームをクローズしてOKです．

        f_in.close();

課題の解答例