前回のあらすじ
熱電流に起因する「雑音」（ノイズ）の除去を行い、夜空の暗い部分を「真っ黒」（つまり階調0）で表現するために、低階調の成分がどの程度の範囲に広がっているか確認した。予想通り背景一面にその成分は広がっていたが、星も10個ほど消えてしまった。果たして何がおきたのか解明する必要がある。そこで、前回はBチャネルのみならず、他のチャネルの様子も見てみよう、というところで終わったのだが、一つやり残したことがあった。それが「ネガとポジ」戦略である。
ネガとポジ
天体観測もかつてはフィルム、そして更にその前は乾板を用いて撮影していた。アンドロメダ銀河中のセファイド型変光星を記録した、有名なハッブルが使用したガラスの写真乾板は「ポジ」(positive)であった。つまり白黒正像画像であった。
今回分析したのは32以下の成分で、あたかも「ネガ」のようであった。あまりにもたくさんの点が該当し問題の分析や処理が大変であった。となれば、その反転画像である「ポジ」で見た方が扱いが簡単になるのは明らかなので、階調32以上の成分をまずは見てみることを思いついた。
やってみたのが、次の画像である。
「ネガ」の画像では「消滅した」と思っていた星像が存在している！しかしノイズらしき点も多数残っていることも確認できた。ということで、もう少しだけバックグランドの値を上げてみよう。試行錯誤の結果、40以上のデータを表示すると「それらしい感じ」となることがわかった。
小さなドットに関しては、どこまでがノイズで、どこからが「暗い恒星」かを区別するのは難しいが、目で認識できる明るめの恒星はすべて表示できていることは（目視で）確認できた。ということで、熱揺らぎに起因するCMOSの「バックグラウンド」は40以下であると仮定することにする。したがって、生データでは255が記録された強度だったとしても、40を引いた215として扱うことにしよう。これは「ダークフレーム処理」に相当する処理である。
それではいよいよ測光である。階調は光の強度の相対値であるから、まずはそれがどんな形になっているか調べておく必要がある。中心部分で明るく、周辺に離れるに従って減光するはずなので、ガウス関数（正規分布）のような形になっていることを期待するわけだが、果たしてどうなっているだろうか？ダークフレーム処理したデータを使って、ミラの中心部分を水平に「切断」した図を見てみたい。具体的にはy=1152の線に沿って、階調をxの関数として表示したものである。
なんとなく、ピークが立っているように見えるのだが、ピークそのものは「ぶった斬られている」ように見える....。その値は215つまり、255-40であり、階調の最高値からバックグランド値を引いたものになっている。これは「飽和」と呼ばれるデジタル画像の症状であることを先ほど知った。「飽和した画像ではピーク部分の光度の情報が失われてしまい、正確な測光ができない」と多くのHPで指摘されている....。たとえば、こちらやこちら。
「飽和」の症状は、特にJPEG形式に一旦変換された画像に出やすいという。Unistellarの場合も、RAWデータは望遠鏡に直接アクセスしてダウンロードしないと手に入らない。通常の方式で画像データをiPadで入手すると、データ送信された段階でpngに変換されてしまい、そこからbmpに戻しても、もはや飽和してしまった情報を取り戻すことはできない。多くのHPで「RAWデータを使おう」とあるが、もはや手遅れである（これからはRAW画像を処理することにする、という方針転換は可能だが、それは未来の話である）。
JPEGというのは、「見た目」が変わらない程度に画像に手を加え圧縮するアルゴリズムである。したがって、明るさが（肉眼では）区別できない「白飛び」の部分は切ってしまおう、という方針なのかもしれない。
しかし、飽和しているといっても根元のデータはしっかり取れているのであるから、この情報を使って、ピーク周辺のデータを復元することは不可能ではあるまい（精度にばらつきはあるだろうが）。形はどうみても正規分布に類似なので、ガウス関数を仮定し、振幅などの係数をフィッティング計算で求めてみた。その結果、次のような復元図（近似形）を得ることができた。
この関数の積分値が「光度」に相当する。ガウス関数に近似したので計算はむしろ楽になった（数値計算が要らない）。 $$ \iint_{\infty^2} A\exp\left(-\frac{x^2+y^2}{d^2}\right)dxdy=A\oint d\theta \int_0^\infty \text{e}^{-\frac{r^2}{d^2}}rdr $$ この面積積分の値は$A\pi d^2$であることは（数物系の人間には）周知の事実である。フィッティング計算の結果は$A=295, d=110$であったから、1月25日のミラの「光度」は（我々の処理においては）
% echo "a=295; d=110; a*3.141592*d^2"|bc -ql 11213912.644000
となった。しかし、物理単位がこの計算にはないので直接他のデータと比較することができない。そのためには、等級が知られている星の光度を上と同じ方法で計算し、その比、つまり相対光度を使って調べる必要がある。次回はその計算をやってみよう。
