Mathematica 4.2でCCDのイメージリダクションを実行するにはどうすればよいのですか．

この例では，Mathematica 4.2がFITS形式で保存されたデータをサポートすることにより，どのような利点があるかを示します．FITSデータには量的なデータが抽出できる天文画像を含んでいることがよくあります．この例で使用するデータは，2000年5月9日にアリゾナ州FlagstaffのLowell観測所にある31インチの望遠鏡を使って収集されたものです．光学システムには，解像度512ピクセルx512ピクセルの液体窒素冷却CCDカメラを使用しました．

重要なシステムの属性はすべて，各画像のヘッダに保存されます．この唯一の例外は，Filter属性の「infrared（赤外線）」が「I」となってしまうことです．ここで使用されたデータは，波長およそ5,500オングストロームで効果のある，半値幅がおよそ890オングストロームである帯域通過のVバンドフィルタを使って収集されました．

リダクションフレームのロード

まず，リダクションの対象となる画像をロードします．フレームは，バイアスフレーム（あるいはゼロフレーム），フラットフィールドフレーム，生画像フレームからなります．この場合，画像はVバンドフィルタを使って撮影されたので，必要なフラットフィールドフレームはVバンドフィルタを通して撮影されたものだけです．すべての場合において，ConversionOptionsはヘッダデータも読み込むよう指定します．

このセクションに含まれるバイアスフレームおよびフラットフィールドフレームは，一連のフレームとは別のセッションで作成されたものです．平均バイアスフレームとは，10枚の異なるバイアスフレームの平均で，平均ラットフィールドフレームとは，Vバンドフィルタで撮影された2つの異なるフラットフィールドフレームの平均です．

バイアスフレームとはCCDチップの"ゼロ"ポイントの測定のためのものです．これは通常露出ゼロで撮影してチップからピクセル値を読み出すことで得られます．リダクションが行われた画像でピクセルカウントを実行する際に，CCD自体によって生成された電子シグナルを含まないようにするためにバイアスフレームが必要となります．バイアス信号の変動はCCDチップ全体でも小さいので，平均バイアスをプロットすると，黒く見えるほとんど一様なフィールドになります．このため，平均バイアスフレームはここではプロットしません．

avgbias = Import["avgbias.fit",
    ConversionOptions -> {"Verbose" -> True}];

次に平均フラットフィールドを読み込んで表示します．正規化されたフラットフィールドフレームは，チップ全体でのCCD感度の不規則性を補うために使用されます．この不規則性は，チップに固有のもの，あるいは，チップの面または使用されているフィルターに付着した塵のためである可能性があります．CCDを一様に光の当たっている表面に露出すると，結果の画像にはこのような不規則性が現れます．このフィールドを正規化して，バイアスを修正した画像を正規化されたフラットフィールドで割ると，画像は補正されます．感度がそれほどよくない部分のピクセルは，チップで最も感度のよい部分のピクセル値と整合性を持つように補正されます．これと同様に，最も感度のよい部分はそれほど感度のよくない部分と整合性を持つように補正されます．

avgvflat = Import["avgvflat.fit",
    ConversionOptions -> {"Verbose" -> True}];

ListDensityPlot[avgvflat[[1,2]],
    Mesh -> False, ImageSize -> 350];

ダークフレームとは，熱で励起された電子によって生成される熱雑音の測定のためのものです．この場合，CCDチップは液体窒素で冷却されるので，感知できるほどの熱カウントはありません．従って，ダークフレームは省略しました．

画像フレームは巨大な楕円銀河M87の画像です．生の画像には右側に幅50ピクセルのオーバースキャン領域（画面からはみ出した部分）が含まれています．この領域は後で削除します．この画像フレームの描画は，次のセクションで説明します．

image = Import["May08_0132.fit",
    ConversionOptions -> {"Verbose" -> True}];

ヘッダを見ることもできます．ヘッダには，画像の種類，使用されたフィルタなどの重要な情報が含まれています．

TableForm[{"OBJECT", "BITPIX", "EXPTIME",
    "DATE-OBS", "OBSERVAT"}/.image[[1, 1]]]

Value -> 'M87'
Value -> 16.
Value -> 300.
Value -> '09/05/00'
Value -> 'lowell'

画像のリダクション

画像のリダクションとは，生画像フレームを分析に利用できるよう準備する，生画像フレーム処理のことです．この例の場合，処理としてすべての画像からのバイアスフレームを差し引き，補正された画像を正規化したフラットフィールドフレームで除算します．

オーバースキャン領域が削除されたらピクセルも削除されるという事実を反映するように元のヘッダを編集するために，元のヘッダは保存しておく必要があります．

oldheader = image[[1, 1]];

オーバースキャン領域を削除するので，ヘッダを編集してピクセル数が減少することを明示しておきます．

newheader = oldheader/.("NAXIS1" -> _) ->
    ("NAXIS1" - > {"Value" -> 512., "Comment" -> ""});

フラットフィールドフレームからオーバースキャン領域を削除し，バイアスフレームも差し引く必要があります．

flatfielddata = Transpose[Take[Transpose[
    avgvflat[[1, 2]]], 512]] - avgbias[[1, 2]];

同様に，生画像フレームからもオーバースキャン領域を削除しなければなりません．

olddata = Transpose[Take[Transpose[
    image[[1, 2]]], 512]];

フラットフィールドフレームを正規化するためには，各ピクセルを全フレームの平均のピクセル値で割ります．

normalizedflatfield = flatfielddata/
    ((Plus @@ Flatten @ flatfielddata)/(512*512.));

ここで，バイアスフレームを差し引いて，正規化したフラットフィールドフレームで割ることにより，画像のリダクション処理を実際に行います．

newdata = (olddata -
    avgbias[[1, 2]])/normalizedflatfield;

最後にデータをヘッダセクションとデータセクションに復元します．このようにしておくと，後で画像を新しいFITS形式にエキスポートしたいときに便利です．

   workingimage = {{newheader, newdata}};

結果の描画

以上で計算はすべて終了したので，結果を見てみましょう．まず，生画像フレームを見ます．画像には"ドーナツ"状のものがたくさんあります．これは，光路ある塵の粒子で，焦点が合わずぼんやりと写っているのです．また，オーバースキャン領域が画像の右側に見えます．オーバースキャン領域には画像データがないため，その部分は前のセクションで削除しました．

ListDensityPlot[image[[1, 2]],
    Mesh -> False, ImageSize -> 350,
    AspectRatio -> Automatic];

これで，リダクションされた画像フレームが表示されます．塵の粒子はフラットフィールディングで削除され，オーバースキャン領域も削除されました．

ListDensityPlot[workingimage[[1, 2]],
    Mesh -> False, ImageSize -> 350];

最後のステップとして，PlotRangeを使うことができます．これは，基本的には画像のコントラストを調整するためです．この場合，M87の中心の超大質量ブラックホールから物質が噴出しているのが分かります．

ListDensityPlot[workingimage[[1, 2]],
    PlotRange -> {75, 2500}, Mesh -> False, ImageSize -> 350];

この例で使用された3つのFITS画像がダウンロードできます．