画像処理ソフト開発ごっこ

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

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

あぷらなーとは基本的に「新しいソフトへの順応性」が低いです。
特に、一度慣れた環境から脱却するのは非常な苦痛を伴います。
ですから天体画像用のソフトは昔からステライメージのみをメインに据えてきていたのですが、そのステライメージ自体、6.5から7への変化についていけなくて、いまだに6.5を多用しています。その6.5ですら各種の謎、たとえば
●ダーク補正時に生じた負の値はどう処理されてる?
●フラット補正時にはなんらかのぼかしフィルタ処理がなされている?
●基準星からのサーチ範囲は、相対的?絶対的?
●割り当てるCPUのコア数は少ない方が高速なケースがある?
●最大エントロピー法で画像が鏡像に変化するケースがある?
●加算比較明コンポジットの具体的なロジックは?
●どうしてデモザイク時に外周部3ピクセルが捨てられる?
●無加工でニコンのNEFファイルをFITSに変換するにはどうする?
などなどを解明するのに何年も掛かる始末でして、なかなか先に進めません。
まあ、その謎解きの過程の中で
『イーブンオッドコンポジット法』とか
『クールファイル補正法』とか
『手動ダーク減算法』とか
『コスミカット法』とか
諸々の新手法を思いついたので、それはそれで楽しいのですが、とにかく個人的には「ブラックボックス化された処理は怖い」んですね・・・・。

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

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


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

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

天体画像処理ソフト自作ごっこ①_f0346040_21271681.jpg
天体画像処理ソフト自作ごっこ①_f0346040_21273445.jpg
数ヶ月にわたる格闘の末(我ながら才能無いなぁ・・・)ついに、FITSファイルの画像データは下記の構造を持っているらしいことが判明しました。

①ヘッダ部分はASCIIコードによる1バイト×80のブロックで1項目(パラメータ)を記述している
②①のデータが合計18セット並んでいる
③ヘッダ部分の直後には0x20(スペース)が②と同じサイズで18セット並んでいる
④次に画像データが格納されている

さらに、④の画像データを読み出すためには
1)最上位ビットが1の場合のみ上位バイトの値から256を引く
2)そこに256を掛ける
3)そこに下位バイトの値を加算する
4)32678を加算する
という演算で1ピクセル分の輝度が復元できることが分かりました。
天体画像処理ソフト自作ごっこ①_f0346040_21390327.jpg
これで、ようやくFITSファイルを表示することができるようになったわけですが・・・・
天体画像処理ソフト自作ごっこ①_f0346040_21400698.jpg
いかんせん素人が力技で演算させているものだから、遅い遅い
なんと1コマのFITSファイルを表示させるのに10~30秒もかかってしまうというゴミアプリが誕生してしまいました。


★奇跡の言語MATLAB降臨

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

だって、ねぇ。
Delphiでこんなに苦しいコード↓を書いてたFITS画像の読み取りが・・・・・
天体画像処理ソフト自作ごっこ①_f0346040_21594996.jpg
MATLABなら

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

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

で終わりですよ?!

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

data_av=(data1+data2)/2;

で終わりですよ?!
驚くべき事に、Forループを回すことすら不要なんです。

誤解を恐れずに言えば、MATLABの変数とは、EXCELで言うところのセルではなくてシートやブック丸ごとなんですね。
これはもう、天体画像処理のために生まれてきたような言語(ちがう)ではないかっ!!
というわけで、画像処理ソフト自作ごっこが一気に加速しました。

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


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

ノータッチガイドを愛するあぷらなーとにとって、コンポジット時の自動位置合わせは非常に大切な処理です。
各種の怪しげなオリジナルロジックを実装する前に、まずはこの自動位置合わせコンポジットをMATLABで実装せねばなりません
もちろん、MATLAB自体は非常に膨大なツール群が用意されていて、課金すればたいていなんでも作れるパーツが手に入ります。ただ、それだと面白くないので、下記のような方針でオリジナルな位置合わせロジックを実装してみることにしました。

①基準星にアンカーを打つ必要すら無いものを目指す
②並進誤差のみを補正すれば充分とする
③できるだけ高速演算を目指す

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

イメージとしては、下記の図のように・・・
天体画像処理ソフト自作ごっこ①_f0346040_22313738.jpg
水色(画像1)と赤(画像2)が重なっていない量を求め、それが最小値を取る状態を「位置合わせできた」と判定するという発想です。上の図で言えば左の状態よりも右の状態の方が位置合わせ精度が高いと判定します。具体的には、画像1と画像2の輝度差分の自乗を積算し、これを位置のズレのパラメータとして使用するという発想ですね。

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

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

★グリッドサーチ法の弱点は「点像」
さきほどの「グリッドサーチ法」は、そもそも連続的な曲面を持つ関数と、その曲面の任意位置から抽出した複数のサンプルデータとをフィッティングするための便法です。したがって、位置合わせする対象(というか輝度分布)が離散的なもの(恒星像などの点像)だと、誤った位置に収束してしまいます。
天体画像処理ソフト自作ごっこ①_f0346040_23065823.jpg
要するに恒星像では位置合わせが不能という意味です。したがって、位置合わせの基準として画面中に存在する星雲などの面積体を自動的にサーチする必要があります。

そこで、まず元画像全体を縦横8×8の64エリアに分割し、それぞれのエリアの輝度メジアンを算出してみました。
天体画像処理ソフト自作ごっこ①_f0346040_23254612.jpg
ここで変数SSは撮像センサーのX方向のピクセル数、ttはY方向のピクセル数です。
たとえば、下記のような元画像なら・・・
天体画像処理ソフト自作ごっこ①_f0346040_23272127.jpg
上記のコードを通すと、こんな感じのエリア判定が得られました。
天体画像処理ソフト自作ごっこ①_f0346040_23281181.jpg
この数値が最も高いエリアを位置合わせ用の基準として判定します。
まずは、採用したエリア以外にマスクを掛けて「のぞき窓」を作るイメージでコードを書きます。

天体画像処理ソフト自作ごっこ①_f0346040_23470942.jpg
やっている演算イメージは、こんな感じです。
天体画像処理ソフト自作ごっこ①_f0346040_23474780.jpg
ちょうど、位置合わせ用の元画像が写ったフィルムをトリミングして、その周囲に黒いイーゼルマスクを掛けた状態ですね。
では、次に位置合わせ対象画像のシフト処理に移ります。


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

まずは、センサーシフト式手ぶれ補正カメラをイメージして、対象画像の周囲に最大シフト幅分の余白を含む平面を定義します。
次に、その平面内を対象画像が丸ごと上下左右にシフトするような仕掛けを準備します。

天体画像処理ソフト自作ごっこ①_f0346040_00130227.jpg
そして、先ほどの周囲をマスキングした基準画像と重ね見るイメージでズレを評価します。
天体画像処理ソフト自作ごっこ①_f0346040_00281435.jpg
さて、次はいよいよ「グリッドサーチ法」の実装です。

★グリッドサーチ法の実装
今回の記事を書いている最中にふと気になってググってみると、近年は機械学習などで用いられているようですね。
ちなみに大昔に書いたあぷらなーとの修士論文を紐解くと、こんな感じで概念が説明されています。
天体画像処理ソフト自作ごっこ①_f0346040_00344310.jpg
要するに本来ならサーチエリア全面に広がるメッシュ状のサーチポイントを、演算速度向上のために9点のグリッドに限定し、徐々にグリッド幅を縮めていくことでベストポジションをサーチするという便法です。
これをMATLABで書くと、こんな感じになりました。
天体画像処理ソフト自作ごっこ①_f0346040_00434472.jpg
ちなみに、1コマの対象画像を基準画像に位置合わせ完了すると、次のコマは前のコマの最適シフト位置を中心として画像シフトが始まるように工夫しています。これにより、多数の画像をコンポジットする際に後半のコマが元画像から大きくズレた場合でもグリッドサーチルーチンが画像をロストするリスクを軽減しています。

★果たしてSI7を越えられたのか??

恐らく、グリッドサーチ法を用いて星雲基準の位置合わせコンポジットを実現した例は希有(皆無?)かと思うのですが、肝心の位置合わせ精度と処理速度はどうでしょうか?
早速、ステライメージ7の「画像マッチングコンポジット」と比較してみましょう。
ASI294MM-Pro+SE120+Hαフィルタで30秒露光したハート星雲の64コマコンポジットで勝負です。



すると・・・
ででん!!
天体画像処理ソフト自作ごっこ①_f0346040_01140767.jpg
※左:SI7でコンポジット 右:自作プログラムでコンポジット

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

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

SI7:Ryzen7 3700X のCPUを100%使い切って200秒
自作:Ryzen7 3700X のCPUを55%使って66秒

ふはははは。
どうやら勝ったな。

MATLABで書いた邪悪な自作画像処理ソフト。
自分にとってはこれこそが正義であることから名付けて
『邪崇帝主』(ジャスティス)じゃぁーッ!!

★★★お約束★★★
①あぷらなーとはプログラミングの素人です。ソースコードの一部が載っていますが、細かいツッコミはご容赦ください。
②お詳しい方は精度の向上や速度向上の策が丸見えでしょうが、「カンニング無しで、試行錯誤する過程が生きがい」という酔狂な性分ゆえ、『正解』は伏せたまま、どうか暖かく見守ってください。
③この続きは不定期に記事にする予定です。

Commented by にゃあ at 2020-12-21 18:35 x
コンポジットのロジックって難しそうだなあと思っていたのですが、なるほど点ではなく面を基準として、そしてその面をさらに狭めていくわけですね。私にはそんな考えがないので、輝点を総当りするのかと思ってましたよ(汗)。しかも、ステライメージより3倍速いとかいうのも愉快でなりません(しかもCPU負荷は半分)。修論がここで生きてくるとか、研究のテーマが何十年もぶれてないのがすごいです。まさに邪崇帝主!
Commented by supernova1987a at 2020-12-21 18:50
> にゃあさん
本来は恒星を基準点にするべきなんですがイマイチ上手く行かなかったので、発想の転換で星雲を基準にして位置合わせしちゃいました。一見、暗くてボヤけてて位置が求まらなさそうに感じますが、恒星と異なり面積がデカい分だけ輝度情報が豊富なので上手く行きました。
邪道バンザイ♪
Commented by taka_ume at 2020-12-24 22:28 x
なるほどこういう方法があるんですね。私の場合、画像を2値化して、OpenCVのテンプレートマッチングを使ったのですが、2値化の閾値は手動で決める仕様なんで、ちょっとめんどくさく。これは是非参考にさせていただきたいです。画像間の平均輝度差も気にしなくて良さそうですし。素晴らしいです。
Commented by supernova1987a at 2020-12-26 17:55
> taka_umeさん
他に例を見ない(?)邪悪なアイディアではありますが、意外と上手く位置合わせできました。これくらい単純な処理ならば、高度なライブラリを用いなくてもサクッとコードが書けちゃいますね♪

※本当は、水面下で色んな手法を次々と試しては玉砕を繰り返していて、最終的にこのアイディアを閃いた・・・・という泥臭い過程がありました。結局、やっていることは、目視で位置合わせする際の「差の絶対値プレビュー」そのものですね。目視で容易に位置合わせできるということはPCにやらせてもいけるのでは?というのが、発端でした。
名前
URL
削除用パスワード
by supernova1987a | 2020-12-21 01:33 | 画像処理ソフト開発ごっこ | Comments(4)

あぷらなーとの写真ブログ


by あぷらなーと