新春から「ででん！！」三連発

あけましておめでとうございます！

今年も邪悪なネタをお届けしますので、よろしくお願いします♪

★今度こそDELPHIと決別の時？

BASIC→C言語→FORTRANと進んだ学生時代。ま、当時はUNIX上でやるシミュレーションや解析の「計算速度が命」でしたので。

社会人になってからはコーディングからは完全に離れていたけれど、ある時にどうしてもWindowsアプリを書かないといけない状況に陥ってしまいDelphi7を購入しました。でも、そのお仕事が終わると、またコーディングを全然やらなくなって１０数年が経過。

天文を再開してしばらくしてZWOの冷却CMOSカメラASI1600MC-Coolに関する諸々の『謎』を解明すべくDelphiの封印を解いてみたものの、とにかくFITSファイルの読み込みが大変で、読み込み用のコードを完成させるだけで３ヶ月かかるという始末。だってエンディアン形式とかすっかり忘れてるし２の補数表現とか未知の領域なのでバグの量産。しかも、素人が書いたコードなので読み込み速度が遅い遅い。１画像読み込むのに30秒もかかるという罰ゲームのようなアプリを作ってしまったという顛末。

それでも、そのコードを用いることでASI1600シリーズについては、色々な発見があったのですが、そろそろ限界。

そこでC＃やPythonに手を出してみたけれど、イマイチしっくりとこない。というわけで、DelphiとExcelを併用した『解析ごっこ』を細々と続けていたのですが・・・・・。

年の瀬に、（全然緩くない『ゆるーく』というパワーワード？で著名な）LambdaさんとTwitter上で議論遊びしていたとき「MATLAB」という聞いたこともない言語名が・・・。調べてみると、かなり昔から理系研究者に愛用されていた解析＆数式処理用の言語らしい。しかも標準でFITSの読み書きに対応しているという素敵さ。ちなみにその特徴は変数が「行列」であることで、ループを回さなくても配列計算ができてしまうとのこと。あの「ほしぞLoveログ」のSamさんも20年ほど使っているそうで、興味が湧いてきたのですが、問題はそのお値段。

Lambdaさん曰く『高級外車が買えるくらい』だそうで、さすがにそれはボツですよねぇ。だって、それだけ軍資金があればアクロマートを100連装しますって（違う）。

ところが、個人用途に絞り込んだ「ホーム」ライセンスってのができてて、そのお値段が・・・・

ええっ！

い・・・１万５千円

・・・ですと！？

思わず後先考えずにポチってしまいました。

★MATLABスゴいよ

さて、爆安とはいえ有料の言語ですから、元を取らないといけません。

早速自作PC２台にインストールして、その言語仕様を探ります。

MATLABは基本的にインタープリタ言語なので、懐かしのN88BASICと同様、プログラムコードを書いて実行することも、コマンドラインから逐次実行することも可能です。驚いたのは、その行列変数なるもので、たとえば

Aという変数に

１　２

３　４

Bという変数に

３　４

５　６

という値が格納されているとすると、

A＋Bとコードを書くだけで

４　６

８　１０

が計算されます。

つまりループで配列要素を１つずつ回すことなく、（まるでExcelの串刺し計算のごとく）配列全体をガバッと演算させられてしまうのですね（数学でいう行列演算も一発で可能）。

・・・ってことは・・・・

これって、まさに天体写真の画像処理＆解析のためにあるような言語じゃないか！！

だって、A,B,Cの各変数に元画像を３枚ロードすれば、

（A+B+C）／３

ってコードを書くだけで加算平均コンポジット完成ですよ！！

これはもう、

「もう、PI行っちゃいな？」

「みんなが使い出したからPIに進みたいけど、真似っこはしない天邪鬼キャラだから進めないのだ」

「やりたいこと、全部できるわよ？」

「いや、そもそも、あぷらなーと家ではフローライトとPIは御法度という家訓なのだ」

などと沼女と押し問答している場合ではありません。

ええ、自分が理想とする画像処理コードをMATLABで書けば良いのですよ！

★格闘すること２日

無い知恵を絞って、比較明コンポジットで生じる光跡の途切れを原理的に軽減する必殺奥義『イーブンオッドコンポジット法』を高速実行するMATLABコードを開発することに成功♪

早速実行してみると、（別ロジックで同等の効果を狙った）ステライメージの「加算比較明コンポジット」で57秒かかるASI1600MCの99コマ処理が、たったの17秒で走る！！

コンパイラではなくインタープリタなのに、スゴいよMATLAB！

★新年の野望を次々撃破

ま、Twitterをご覧になっている方はご存じだと思うんですが、MATLABで下記のような（数年間温めてきた）プロジェクトが一気に実現しそうです。

①ダークを引くと背景は荒れるのか？の検証

例えば、ダークを120コマ確保していたとして、それを色んな枚数でコンポジットして複数のダークファイルの「効き」を比較するのは難儀なものです。そこで、色んな条件のダークファイルを一発で自動作成するコードをMATLABで書いてみました。

次に、その色んなダークファイルを用いてSI6.5で仕上げた画像を一気に読み込んで、星が写っていない背景を切り取ってノイズ量（正確には輝度の空間揺らぎの標準偏差）を測定するMATLABコードを書きました。

すると・・・・・

ででん！！

※左：ダーク無し　中：ダーク１コマ適用　右：ダーク120コマ適用

まず、このように、ダークを１枚適用すると（無論ホットピクセルは消えるものの）背景の縮緬ノイズはかえって悪化するが、十分な枚数のダークを引いた場合は、縮緬ノイズも消滅することが確かめられました。

　　　※対数グラフで演算値が負に発散しないように、ダークを作用させないケースを「ダーク枚数0.1枚」として扱ってあります

ただし、上のグラフをみると分かるように、これまでの経験値よりも少ない枚数でその効果が飽和してしまっているように見えます。ところが、実際にそれらの画像を目視で比較してみると

※左：ダーク16枚適用　右：ダーク120枚適用

このように、明瞭な「効き」の差があります。つまり、一般的な「画像内の輝度揺らぎの標準偏差を測定する」方法では縮緬ノイズは評価できないことになります。

ここから、以前から主張している怪しい理論「同じノイズ量を持つ画像でも『心理的エントロピー』によって、それを人間がノイズとして認識できるかどうかが決まる」説が有力になりました。

②『心理的エントロピー』を数値化する糸口

ちょうど、次のコードを書いていたときにLambdaさんから「心理的エントロピーとは、結局、空間周波数の事なのではないか？」とのご指摘を受けたのですが、実は私もうすうす感じていてFFT（高速フーリエ変換）の実装に取りかかっていたのです。

ところが、これ（周波数成分）が、うまく検出できません。

その後、色々とアドバイス（FFTにかける前に全体の平均輝度を減じればどうか）を頂いたので、それらを実装してみました。

すると・・・・

ででん！！

明瞭に高周波成分（縮緬ノイズが走る方向と直交して生じるシマシマ）を検出することに成功しました。

次に解析の妥当性を検討するため、同じ画像のうち、恒星が複数写っている領域を位置合わせ無しでコンポジットし、同様のFFTを掛けてみると

このように、同様の周波数成分が検出されました。

あとは、これを数値化するロジックを構築すれば、縮緬ノイズの機械的な評価パラメータが見出せるかも知れませんし、『手動ダーク減算法』や『クールファイル補正法』でも消しきれなかった縮緬ノイズ軽減の糸口が掴めるかもしれません。

③不良ピクセルを排除する『ソフトウェア・ピクセルマッピング法』の開発糸口

以前、Delphiやマカリで画像解析していた際、どうやらカメラ内のピクセルは、その個々によって特性にバラツキがあり、平均輝度値がズレているピクセル（これは有名で、ダーク減算やバイアス減算で軽減しますよね）の他に、輝度値が酔っ払いのようにふらつく『酩酊ピクセル』が混在していそうだ、という仮説を立てていました。

　　　※上：通常ピクセル　下：『酩酊ピクセル』

このように、平均輝度が同じなのに輝度変化（時間的ゆらぎ）が全く異なるピクセルが混在しているのですね。

この『酩酊ピクセル』が存在していると、画像処理後に残存するノイズの原因になるばかりか、重要なダークファイルやフラットファイルの有効性を破壊してしまう恐れがあります。この不良ピクセルをソフトウェア的にカットするのが『ソフトウェア・ピクセルマッピング法』です。そのためには、カメラ内の全ピクセルについて、その平均輝度値とその揺らぎを測定する必要があります。

さすがに、これは手作業では不可能な膨大なお仕事です。

そこで、MATLABを用いて、「３２コマのダークファイルを読み込み、３２層の３次元FITSとして格納」し、それを個々のピクセル毎に解析することによって、「各ピクセルの性質を全数調査する」というコードを書いてみました。

ASI1600MMは１ファイルで約１６３９万個のピクセル。それが３２層です。つまり「３２個のサンプルデータから平均値と標準偏差を解析する」調査を「１６３９万回」こなす・・・・・こんな邪悪な作業量なのに、MATLABはたったの５７秒でその全てを演算してくれました。

すると・・・

ででん！！

ふはははは。

よもや、こんな邪悪な解析を行うヒマなアマチュア天文家はおるまい。

これは、非常に面白い事になってきました。

平均輝度と標準偏差のスキャッタプロットを見ると明らかなように、どうも想像以上に色んな特性のピクセルが混在しているようです。

素人なので、勘違いしている可能性も高いですが、例えば、上のように特性が異なる５～６群のピクセルに分類できそうです。

特に画像処理をする上で致命的な結果を招きかねないのが『酩酊ピクセル』（グラフ最上部の群）でして、対数グラフで見てるとそうでもないですが、これをリニアスケールで３Dプロットしてみると

※輝度値では無く揺らぎ量（標準偏差）のプロットである点にご注意

各画素の輝度揺らぎ（輝度標準偏差）は平均して257でした。それに対して、上のグラフでは実にその10倍以上の揺らぎを持つ『酩酊ピクセル』がウジャウジャ存在していることが分かります。ちなみに、平均して輝度がズレるピクセル（ホットピクセルやなど）ならば、ダーク減算などで補正することが可能です。しかし、輝度揺らぎ（一種の時間ノイズ）は原理上コンポジット枚数を増やすしか軽減方法がありません。今回検出した『酩酊ピクセル』は平均的なピクセルと比べその輝度揺らぎが10倍以上もあるのですから、理論的にはコンポジット枚数を100倍（10の2乗）に増やさないと同等の結果が得られません。逆に言えば、『酩酊ピクセル』はせっかく稼いだコンポジット枚数を100分の１に減らしてしまう凶悪なピクセルだと言えるでしょう。

さて、今回の解析で、重要なのはこれらの解析データが素子全体の「グループ」としてではく、個々のピクセル毎に保存されている点です。したがって、今後バイアスの特性やゲインの特性や温度の特性などの検査をすることによって、どのピクセルが『悪さ』をしているかその『犯人』が特定できるということです。画像処理の方向性には色々な正解があるのですが、あぷらなーとは「できるかぎりのキャリブレーションを行う」ことで「使い物にならないピクセルは排除」してから画像処理に回し、「必要最小限のレタッチに止める」ことを目指したいと思います。

そのための秘策が『ソフトウェア・ピクセルマッピング法』というわけです。

★というわけで

令和２年一発目のエントリーでは「新年の邪悪な抱負を列挙」するだけのつもりだったのに、MATLAB降臨によって、3日間で１年分くらいプロジェクトが進んでしまったよーというお話でした。

この邪悪な記事を見てゾクゾクされた方(いるのか？）は・・・・

「一家に一台・MATLAB」！

・・・かもね♪

※追記

対象画像ファイルを120コマに変更して解析した最新のグラフを追加します。

各群の特徴を分かりやすくするため、片対数グラフで貼り付けておきます。

※グラフのキャプションに32コマとあるのは120コマの誤記です。

※追記その2

『酩酊ピクセル』についてさらに解析を進めると、法外な標準偏差でランダムに変動しているのではなく、突発的に異常値を吐いた結果であることが判明。同様の傾向はASI294やASI174でも見られるため、特性ではなく宇宙線のヒットなどの現象ではないかと推測。トリウムレンズの放射線に被曝させて、現象をある程度再現できたかのように見えました。ところが、さらに色々なデータを見てみると、宇宙線のヒットだとするにはその分布があまりにもタイトである点が気にかかり、現在さらなる解析を実行中です。もう少し時間を下さい。