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

★新しいソフトへの順応性は低いです

あぷらなーとは基本的に「新しいソフトへの順応性」が低いです。

特に、一度慣れた環境から脱却するのは非常な苦痛を伴います。

ですから天体画像用のソフトは昔からステライメージのみをメインに据えてきていたのですが、そのステライメージ自体、6.5から７への変化についていけなくて、いまだに6.5を多用しています。その6.5ですら各種の謎、たとえば

●ダーク補正時に生じた負の値はどう処理されてる？

●フラット補正時にはなんらかのぼかしフィルタ処理がなされている？

●基準星からのサーチ範囲は、相対的？絶対的？

●割り当てるCPUのコア数は少ない方が高速なケースがある？

●最大エントロピー法で画像が鏡像に変化するケースがある？

●加算比較明コンポジットの具体的なロジックは？

●どうしてデモザイク時に外周部３ピクセルが捨てられる？

●無加工でニコンのNEFファイルをFITSに変換するにはどうする？

などなどを解明するのに何年も掛かる始末でして、なかなか先に進めません。

まあ、その謎解きの過程の中で

『イーブンオッドコンポジット法』とか

『クールファイル補正法』とか

『手動ダーク減算法』とか

『コスミカット法』とか

諸々の新手法を思いついたので、それはそれで楽しいのですが、とにかく個人的には「ブラックボックス化された処理は怖い」んですね・・・・。

★ならば、全部自分で作っちゃえ

というわけで、いっそのこと「自分に必要な機能」だけを自力でプログラミングしちゃえ！と無謀なプロジェクトがスタートしました。かれこれ４年ほど前のことです。

★まずはFITSの読込ができないとお話にならない

なにはともあれ、まずはFITSファイルの読み込みという難関をクリアせねばなりません。

FITSファイルに対応していないDelphiでコードを書いていた頃は、FITSファイルのデータを１バイトずつ読み込んで、どこに何が格納されているのかを手探りで推測するという荒行を行っていました。

数ヶ月にわたる格闘の末（我ながら才能無いなぁ・・・）ついに、FITSファイルの画像データは下記の構造を持っているらしいことが判明しました。

①ヘッダ部分はASCIIコードによる1バイト×８０のブロックで１項目（パラメータ）を記述している

②①のデータが合計１８セット並んでいる

③ヘッダ部分の直後には0x20(スペース）が②と同じサイズで18セット並んでいる

④次に画像データが格納されている

さらに、④の画像データを読み出すためには

１）最上位ビットが１の場合のみ上位バイトの値から２５６を引く

２）そこに２５６を掛ける

３）そこに下位バイトの値を加算する

４）３２６７８を加算する

という演算で１ピクセル分の輝度が復元できることが分かりました。

これで、ようやくFITSファイルを表示することができるようになったわけですが・・・・

いかんせん素人が力技で演算させているものだから、遅い遅い。

なんと１コマのFITSファイルを表示させるのに１０～３０秒もかかってしまうというゴミアプリが誕生してしまいました。

★奇跡の言語MATLAB降臨

一時は頓挫しかけたこのプロジェクトですが、ちょうど１年ほど前にTwitter上の会話からLambdaさんにMATLABという言語を薦めていただいて触れてみたところ・・・・・「なんじゃ、こりゃ？！」というほど衝撃的な言語仕様。まさに「こんな言語が欲しかった！」とヒザを打ちました。

だって、ねぇ。

Delphiでこんなに苦しいコード↓を書いてたFITS画像の読み取りが・・・・・

MATLABなら

data = cast(fitsread(myfilename),'int16');

※myfilenameは対象画像のフルパス

で終わりですよ？！

しかも、data1という変数名の画像（１枚の画像が丸ごと格納されている2次元行列変数）とdata2という変数名の画像を加算平均コンポジットするコードなんて、

data_av=(data1+data2)/2;

で終わりですよ？！

驚くべき事に、Forループを回すことすら不要なんです。

誤解を恐れずに言えば、MATLABの変数とは、EXCELで言うところのセルではなくてシートやブック丸ごとなんですね。

これはもう、天体画像処理のために生まれてきたような言語（ちがう）ではないかっ！！

というわけで、画像処理ソフト自作ごっこが一気に加速しました。

※この辺の顛末は、恐れ多くもMATLAB公式さんのサイトに掲載↓していただいてます♪

★最初の難関は「位置合わせコンポジット」

ノータッチガイドを愛するあぷらなーとにとって、コンポジット時の自動位置合わせは非常に大切な処理です。

各種の怪しげなオリジナルロジックを実装する前に、まずはこの自動位置合わせコンポジットをMATLABで実装せねばなりません。

もちろん、MATLAB自体は非常に膨大なツール群が用意されていて、課金すればたいていなんでも作れるパーツが手に入ります。ただ、それだと面白くないので、下記のような方針でオリジナルな位置合わせロジックを実装してみることにしました。

①基準星にアンカーを打つ必要すら無いものを目指す

②並進誤差のみを補正すれば充分とする

③できるだけ高速演算を目指す

例えるなら、ステライメージに実装されている「自動画像マッチング」の超高速版ですね。

イメージとしては、下記の図のように・・・

水色（画像１）と赤（画像２）が重なっていない量を求め、それが最小値を取る状態を「位置合わせできた」と判定するという発想です。上の図で言えば左の状態よりも右の状態の方が位置合わせ精度が高いと判定します。具体的には、画像１と画像２の輝度差分の自乗を積算し、これを位置のズレのパラメータとして使用するという発想ですね。

恐らく、一般的にはメッシュ法（一方の画像を必要とされる位置合わせ精度と一致する間隔を持つ格子点上に逐次配置していき、自乗誤差を評価して最小値を探る）か、任意方向について評価誤差量変化のグラディエント（微分係数）を求めていき、それがゼロになった所を極小値であると判定する方法が用いられているのではないかと想像するのですが、ガイドエラーが大きい場合は膨大な演算時間が必要そうです。

こういうときは、やはり「昔取った杵柄」の投入でしょう（笑）。

学生時代に取り組んでいた「高エネルギー１次宇宙線が地球大気と相互作用した結果生じる２次宇宙線粒子群」が地表に到達した際の「粒子数密度のラテラル分布から入射位置を高速推定」するための必殺技「グリッドサーチ法」を天体画像用にアレンジすることに挑戦です。さらに、その際の位置合わせ作業としては「２コマのポジフィルム像を手動で重ねる必要が生じたら、どういう方法をとるか」をイメージしました。

★グリッドサーチ法の弱点は「点像」

さきほどの「グリッドサーチ法」は、そもそも連続的な曲面を持つ関数と、その曲面の任意位置から抽出した複数のサンプルデータとをフィッティングするための便法です。したがって、位置合わせする対象（というか輝度分布）が離散的なもの（恒星像などの点像）だと、誤った位置に収束してしまいます。

要するに恒星像では位置合わせが不能という意味です。したがって、位置合わせの基準として画面中に存在する星雲などの面積体を自動的にサーチする必要があります。

そこで、まず元画像全体を縦横８×８の６４エリアに分割し、それぞれのエリアの輝度メジアンを算出してみました。

ここで変数SSは撮像センサーのX方向のピクセル数、ttはY方向のピクセル数です。

たとえば、下記のような元画像なら・・・

上記のコードを通すと、こんな感じのエリア判定が得られました。

この数値が最も高いエリアを位置合わせ用の基準として判定します。

まずは、採用したエリア以外にマスクを掛けて「のぞき窓」を作るイメージでコードを書きます。

やっている演算イメージは、こんな感じです。

ちょうど、位置合わせ用の元画像が写ったフィルムをトリミングして、その周囲に黒いイーゼルマスクを掛けた状態ですね。

では、次に位置合わせ対象画像のシフト処理に移ります。

★位置合わせ対象画像の処理

イメージとしては、上記のような位置合わせ基準画像を中央に貼り付けたイーゼルに対象画像を重ねて、像のズレが最小になるまで上下左右にシフトさせる感じでコードを書いてみます。

まずは、センサーシフト式手ぶれ補正カメラをイメージして、対象画像の周囲に最大シフト幅分の余白を含む平面を定義します。

次に、その平面内を対象画像が丸ごと上下左右にシフトするような仕掛けを準備します。

そして、先ほどの周囲をマスキングした基準画像と重ね見るイメージでズレを評価します。

さて、次はいよいよ「グリッドサーチ法」の実装です。

★グリッドサーチ法の実装

今回の記事を書いている最中にふと気になってググってみると、近年は機械学習などで用いられているようですね。

ちなみに大昔に書いたあぷらなーとの修士論文を紐解くと、こんな感じで概念が説明されています。

要するに本来ならサーチエリア全面に広がるメッシュ状のサーチポイントを、演算速度向上のために９点のグリッドに限定し、徐々にグリッド幅を縮めていくことでベストポジションをサーチするという便法です。

これをMATLABで書くと、こんな感じになりました。

ちなみに、１コマの対象画像を基準画像に位置合わせ完了すると、次のコマは前のコマの最適シフト位置を中心として画像シフトが始まるように工夫しています。これにより、多数の画像をコンポジットする際に後半のコマが元画像から大きくズレた場合でもグリッドサーチルーチンが画像をロストするリスクを軽減しています。

★果たしてSI７を越えられたのか？？

恐らく、グリッドサーチ法を用いて星雲基準の位置合わせコンポジットを実現した例は希有（皆無？）かと思うのですが、肝心の位置合わせ精度と処理速度はどうでしょうか？

早速、ステライメージ７の「画像マッチングコンポジット」と比較してみましょう。

ASI294MM-Pro＋SE120＋Hαフィルタで３０秒露光したハート星雲の６４コマコンポジットで勝負です。

すると・・・

ででん！！

※左：SI7でコンポジット　右：自作プログラムでコンポジット

あんな「冗談だろ？」というほど単純で大胆なロジックなのに、互角の精度を示しているようです！！

そして、その処理速度は・・・・・

SI7：Ryzen7 3700X のCPUを100％使い切って200秒

自作：Ryzen7 3700X のCPUを55％使って66秒

ふはははは。

どうやら勝ったな。

MATLABで書いた邪悪な自作画像処理ソフト。

自分にとってはこれこそが正義であることから名付けて

『邪崇帝主』（ジャスティス）じゃぁーッ！！

★★★お約束★★★

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

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

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