天体画像処理ソフト自作ごっこ②

★新年明けましておめでとうございます

旧年中はお世話になりました。

本年もよろしくお願いします

というわけで、本年の抱負を・・・・

①自作画像処理ソフト『邪崇帝主（ジャスティス）』を完成させる

②自作オンアキシスガイダー『テレサ・改』を完成させる

③『光害チョッパー』を完成させる

④「スペクトロヘリオグラム」を作る

⑤『メカニカルPEC』で星空を実写する

むう。我ながら邪悪すぎるぞ♪

・・・というわけで、まずは①の続きから開始です。

★グリッドサーチ法の活用可能性

前回記事では、グリッドサーチ法による位置合わせ（アライメント）機能の実装について書きました。

ちなみに、年末に発売されたステライメージ９は動作がかなり高速化されており、特にメトカーフバッチを利用した恒星基準位置合わせコンポジットはメモリ消費量も最小限に止められる上に『爆速』であることを確認しました。

こうなると、わざわざ位置合わせコンポジットのルーチンを実装する意義は薄れるのですが、カラー画像をデモザイク（ディベイヤー）前にR/G1/G2/Bに分離して画像処理する手法には、それなりのメリットがありますし、基準星にアンカーを打たない位置合わせ方は画面中に輝点がない面積体（月面とか惑星とか）のスタックにも流用できるというメリットもありますので、開発ごっこは継続です。

使用言語であるMATLAB home　は年末に2019ｂから2020bへバージョンアップさせました。今のところ特に不具合もなく絶好調で動いてくれています。

実際はより高度な処理を行う各種パーツが用意されているのですが、結果重視ではなく過程重視の遊びなので、当面は自力でロジックを考えて泥臭く攻めます。

さて、恒星ではなく星雲の高輝度エリアをズレ判定の基準としたグリッドサーチ法では原理上惑星などのスタックにも転用できそうなので試してみました。これはステライメージに実装されている「画像マッチング」にも言えることです。

VMC260L＋ASI1600MM-Coolを640x480ピクセルでROI（クロップ）して撮像したSer動画ファイルから抜き出した256コマの画像を素材にして

①SI９で画像マッチング

②邪崇帝主でグリッドサーチコンポジット

③AS3でスタック

を試してみましょう。

なお、比較のためAS3の良像選択はOFFにしています。

　　　※左：SI9　中：邪崇帝主　右：AS3

とりあえずグリッドサーチ法でも木星の位置合わせスタックは成功するようです。

次に、それぞれのスタック画像をレジスタックスに掛けてウェーブレット処理してみましょう。

　　　※左：SI9　中：邪崇帝主　右：AS3

良像選択を施さない限りは、三者の間に明確な差は生じないようですね。

さて、処理速度はどうでしょうか。

Ryzen7 3700X+メモリ96GBのPCで256コマのスタック処理に要した時間は

　SI9：87秒

　邪崇帝主：14秒

　AS3：４秒

でした。

さすがに専用ソフトであるAS3の爆速っぷりは別格ですが、自作のロジックもなかなか健闘してると思います。

★加算コンポジットの実装

ここからは、ほとんど自分のための覚え書きになりますが、自作『邪崇帝主』における各種コンポジットの演算ロジックをまとめておこうと思います。もちろん、位置合わせが済んでいるか位置合わせが不要な画像に対するお話です。

まずは、加算コンポジットです。

これは単純で下記のように各コマの同一座標に位置するピクセルの輝度を足し算するだけです。

ただし、MATALABを用いて記述する場合、実際のコードはかなり特殊です。

普通はFORループで１ピクセルずつ読み出して全コマ分を加算していくコードになるのですが、MATALABは変数自体がデフォルトで行列（２次元配列）になっているため、画像内で座標をループさせる必要ありません。

たとえば、

C＝A+B；

とコードを書くだけで、画像A丸ごとと画像B丸ごとの加算コンポジット画像がCに格納されます。後は、これを全コマ数分だけ繰り返すだけなので楽勝です。

★加算平均コンポジットの実装

もっとも一般的なコンポジットである加算平均コンポジットについても単純です。

先ほどの加算コンポジットの結果をコマ数で割り算するだけで終わりです。

Cを加算コンポジット後の画像を格納した配列名、ｎをコマ数とすれば、

D＝C/n;

というコードを書けば終わりです。

なお、手作業でコンポジットする際は、いわゆる『トーナメント方式』と呼ばれる２コマずつコンポジットさせていく方法が必要でした。しかし恐らくそれは先に合計を取ってから枚数で除算しようとすると合計を演算している途中で輝度値がオーバーフローするからだと思います。ただしステライメージやMATALBは内部演算の器が大きい（SIは96bit、MATLABは64bit）のため特にこれを心配する必要はありません。

たとえば、MATALABで16bit出力カメラの画像を加算コンポジットする場合なら

２＾（６４－１６）＝２．８１×１０＾１４

となりますので、何コマでもサチらずにコンポジット可能ですね。

★比較明コンポジットの実装

比較明コンポジットのロジックは対象画像同士の同一座標ピクセルを比較して、そのうち明るい方の輝度を採用するというものです。

最初に思いつくのは全コマをロードしておき、任意のピクセルについて全コマ中の最大値を求めて採用するというものですが、これだと演算中の消費メモリが大きいため１コマずつ演算してメモリを解放する方法を考えてみました。

このロジックのキモはどちらのピクセルが明るいかの判定結果（正なら１、偽なら０）を格納したマスク的な配列を用意し、それを元画像に掛けることで必要なピクセル以外が消えるようにしたことです。

これをMATALABで書くと

こんな感じになりました。

注意すべきは途中に出てくる掛け算です。

MATLABでは、その仕様上

C=A＊B；

と書くとAとBの行列積になってしまうので、これを各要素毎の積にするために

C＝A．＊B；

というように変数名の後にドットを付けることぐらいでしょうか。

★コスミカット法の実装

『コスミカット法』は2次宇宙線などの高エネルギー荷電粒子がカメラのセンサーをヒットした際に生じる突発的なノイズを高速に除去するために考案した手法です。

放射線ヒットノイズは通常の輝点ノイズとは異なり、入射角度が浅いとひっかき傷のような構造を持ちます。したがって、通常のホット除去フィルタ処理では排除が困難です。そこで一般的にはシグマクリップを用いて時系列上の輝度変化の外れ値として排除する方法が考えられます。ただしシグマクリップはメモリの消費も計算コストも高いので、その処理に結構時間がかかり大変です。

そこで、放射線ヒットノイズが他のノイズよりも極めて高い値を持つ（弱い物もあるが悪影響が無い）ことに注目し、時系列上の最大輝度値を弾くという方向性でロジックを組み立ててみました。

上記のAとB中の赤いセルが放射線ヒットが載ったピクセルを示します。

全体的な方向性は極めて単純で、加算コンポジット画像から比較明コンポジット画像を引いて、全コマ数ー１　で割り算するだけです。

この手法の面白いところは無理に自力でコードを書かなくてもステライメージなどの既存ソフトで容易に実現（加算コンポジットと比較明コンポジットの画像を保存しておき、それらを減算コンポジットするだけ）できるところです。ただ、せっかくなので『邪崇帝主』に実装してみましょう。

すでに加算コンポジットと比較明コンポジットのルーチン（水色破線）は実装済なので話は簡単。加算コンポジット画像をA、比較明コンポジット画像をB、撮影コマ数をｎとするとき、

C＝（A-B）／（n-1）；

という主旨のコードをたった1行加筆するだけで完成しました（笑）。

非常に単純な仕掛けですが、結構効果があります。　※左：加算平均　右：コスミカット

上記写真の左は通常の加算平均コンポジットを行った結果で、黄色の丸で囲んだ部分に放射線ヒットノイズが写っています。これが、コスミカット法を用いることで消えます。また、上記のように人工衛星の光跡を除去するのにも効果があります。

★イーブンオッド法の実装

『イーブンオッド法』は、いわゆる星グル写真（固定撮影で星の軌跡を写す）などを撮影する際に比較明コンポジットでは原理的に生じてしまう光跡の途切れを（ボカシなどで誤魔化さずに）解消するために考案したロジックです。

上記のまとめ記事中にも述べたように、このロジックはFlatAideProに実装していただいて表の世界にデビューすることができました。また、全く別ロジックで数学的には全く同様の結果を得る「加算比較明コンポジット」がステライメージには実装されているので、今さら感が強いのですが、思い入れのある手法なので『邪崇帝主』に実装を試みてみましょう。ザックリとしたロジックは下記の通りです。

これをMATLABで書くと

こんな感じになりました。

コマNoを２で割ったときの余りが１か０かでオッド群とイーブン群に弁別して、上記の水色破線のルーチンでそれぞれ比較明コンポジットします。最後にそれらを赤破線の行で加算しておしまいです。

こんな単純なロジックですが、結構効果があります。

実際に、ニコンD610＋シグマ20mmF1.8で固定撮影した15秒露光のRAW画像100コマを処理してみると下記のようになりました。

　　※左：比較明　右：イーブンオッド

このように比較明で途切れていた光跡がイーブンオッドだとつながります。

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

今回は、ここまで。

こんな「だれ得」な記事にお付き合いいただいた方、ありがとうございました。

今後もちょくちょく続きを書きます。

★★★お約束★★★

①あぷらなーとはプログラミングの素人です。ソースコードの一部を載っていますが、細かいツッコミはご容赦ください。

②お詳しい方は精度の向上や速度向上の策が丸見えでしょうが、「カンニング無しで試行錯誤する過程が生きがい」という酔狂な性分ゆえ、『正解』は伏せたまま、どうか暖かく見守ってください。

③この続きは不定期に記事にする予定です。