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

あけましておめでとうございます!
今年も邪悪なネタをお届けしますので、よろしくお願いします♪

★今度こそDELPHIと決別の時?
BASIC→C言語→FORTRANと進んだ学生時代。ま、当時はUNIX上でやるシミュレーションや解析の「計算速度が命」でしたので。
社会人になってからはコーディングからは完全に離れていたけれど、ある時にどうしてもWindowsアプリを書かないといけない状況に陥ってしまいDelphi7を購入しました。でも、そのお仕事が終わると、またコーディングを全然やらなくなって10数年が経過。

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

それでも、そのコードを用いることでASI1600シリーズについては、色々な発見があったのですが、そろそろ限界
そこでC#やPythonに手を出してみたけれど、イマイチしっくりとこないというわけで、DelphiとExcelを併用した『解析ごっこ』を細々と続けていたのですが・・・・・。

年の瀬に、(全然緩くない『ゆるーく』というパワーワード?で著名な)LambdaさんとTwitter上で議論遊びしていたとき「MATLAB」という聞いたこともない言語名が・・・。調べてみると、かなり昔から理系研究者に愛用されていた解析&数式処理用の言語らしい。しかも標準でFITSの読み書きに対応しているという素敵さ。ちなみにその特徴は変数が「行列」であることで、ループを回さなくても配列計算ができてしまうとのこと。あの「ほしぞLoveログ」のSamさんも20年ほど使っているそうで、興味が湧いてきたのですが、問題はそのお値段
Lambdaさん曰く『高級外車が買えるくらい』だそうで、さすがにそれはボツですよねぇ。だって、それだけ軍資金があればアクロマートを100連装しますって(違う)。

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

ええっ!

い・・・1万5千円

・・・ですと!?

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


★MATLABスゴいよ

さて、爆安とはいえ有料の言語ですから、元を取らないといけません。
早速自作PC2台にインストールして、その言語仕様を探ります

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

Aという変数に

1 2
3 4

Bという変数に

3 4
5 6

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

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

4 6
8 10

が計算されます。

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

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

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

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

(A+B+C)/3

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

これはもう、
「もう、PI行っちゃいな?」
「みんなが使い出したからPIに進みたいけど、真似っこはしない天邪鬼キャラだから進めないのだ」
「やりたいこと、全部できるわよ?」
「いや、そもそも、あぷらなーと家ではフローライトとPIは御法度という家訓なのだ」
などと沼女と押し問答している場合ではありません

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


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

新春から「ででん!!」三連発_f0346040_05552633.jpg
早速実行してみると、(別ロジックで同等の効果を狙った)ステライメージの「加算比較明コンポジット」で57秒かかるASI1600MCの99コマ処理が、たったの17秒で走る!!
新春から「ででん!!」三連発_f0346040_05514197.jpg

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

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


①ダークを引くと背景は荒れるのか?の検証
例えば、ダークを120コマ確保していたとして、それを色んな枚数でコンポジットして複数のダークファイルの「効き」を比較するのは難儀なものです。そこで、色んな条件のダークファイルを一発で自動作成するコードをMATLABで書いてみました
次に、その色んなダークファイルを用いてSI6.5で仕上げた画像を一気に読み込んで、星が写っていない背景を切り取ってノイズ量(正確には輝度の空間揺らぎの標準偏差)を測定するMATLABコードを書きました

すると・・・・・

ででん!!
新春から「ででん!!」三連発_f0346040_08250430.jpg

※左:ダーク無し 中:ダーク1コマ適用 右:ダーク120コマ適用


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

新春から「ででん!!」三連発_f0346040_06242381.jpg
   ※対数グラフで演算値が負に発散しないように、ダークを作用させないケースを「ダーク枚数0.1枚」として扱ってあります

ただし、上のグラフをみると分かるように、これまでの経験値よりも少ない枚数でその効果が飽和してしまっているように見えます。ところが、実際にそれらの画像を目視で比較してみると
新春から「ででん!!」三連発_f0346040_08292527.jpg
※左:ダーク16枚適用 右:ダーク120枚適用

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

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



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

ところが、これ(周波数成分)が、うまく検出できません。
新春から「ででん!!」三連発_f0346040_06420279.jpg

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

すると・・・・

ででん!!
新春から「ででん!!」三連発_f0346040_06440436.jpg
明瞭に高周波成分(縮緬ノイズが走る方向と直交して生じるシマシマ)を検出することに成功しました。

次に解析の妥当性を検討するため、同じ画像のうち、恒星が複数写っている領域を位置合わせ無しでコンポジットし、同様のFFTを掛けてみる
新春から「ででん!!」三連発_f0346040_06482472.jpg
このように、同様の周波数成分が検出されました。
あとは、これを数値化するロジックを構築すれば、縮緬ノイズの機械的な評価パラメータが見出せるかも知れませんし、『手動ダーク減算法』や『クールファイル補正法』でも消しきれなかった縮緬ノイズ軽減の糸口が掴めるかもしれません。


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

以前、Delphiやマカリで画像解析していた際、どうやらカメラ内のピクセルは、その個々によって特性にバラツキがあり、平均輝度値がズレているピクセル(これは有名で、ダーク減算やバイアス減算で軽減しますよね)の他に、輝度値が酔っ払いのようにふらつく『酩酊ピクセル』が混在していそうだ、という仮説を立てていました。
新春から「ででん!!」三連発_f0346040_07093472.jpg
   ※上:通常ピクセル 下:『酩酊ピクセル』

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

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

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

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

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

すると・・・

ででん!!
新春から「ででん!!」三連発_f0346040_07244152.jpg
ふはははは。

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

これは、非常に面白い事になってきました。
平均輝度と標準偏差のスキャッタプロットを見ると明らかなように、どうも想像以上に色んな特性のピクセルが混在しているようです。
新春から「ででん!!」三連発_f0346040_07272773.jpg
素人なので、勘違いしている可能性も高いですが、例えば、上のように特性が異なる5~6群のピクセルに分類できそうです。
特に画像処理をする上で致命的な結果を招きかねないのが『酩酊ピクセル』(グラフ最上部の群)でして、対数グラフで見てるとそうでもないですが、これをリニアスケールで3Dプロットしてみると
新春から「ででん!!」三連発_f0346040_07412919.jpg
※輝度値では無く揺らぎ量(標準偏差)のプロットである点にご注意

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

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

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


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

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

「一家に一台・MATLAB」!

・・・かもね♪



※追記
対象画像ファイルを120コマに変更して解析した最新のグラフを追加します。
各群の特徴を分かりやすくするため、片対数グラフで貼り付けておきます。
※グラフのキャプションに32コマとあるのは120コマの誤記です。
新春から「ででん!!」三連発_f0346040_09422746.jpg
※追記その2
『酩酊ピクセル』についてさらに解析を進めると、法外な標準偏差でランダムに変動しているのではなく、突発的に異常値を吐いた結果であることが判明。同様の傾向はASI294やASI174でも見られるため、特性ではなく宇宙線のヒットなどの現象ではないかと推測。トリウムレンズの放射線に被曝させて、現象をある程度再現できたかのように見えました。ところが、さらに色々なデータを見てみると、宇宙線のヒットだとするにはその分布があまりにもタイトである点が気にかかり、現在さらなる解析を実行中です。もう少し時間を下さい。


Commented by te kure at 2020-01-06 09:57 x
明けましておめでとうございます🎊
そうでしたか・・ 当初の一見お気楽そうな“ででん!”解析には大変な努力が払われてたのですね〜
超強力な武器 MATLAB は素晴らしいお年玉ですね。
せっかくのお正月🎍休みにご苦労様でした。(笑)

しかし当方 Basic も触った事がありません。
その成果の恩恵に預かれるのやら・・?😅
Commented by kem2017 at 2020-01-06 13:02
こりゃまた難読なブログだこと (^o^;

A= dataEv>=data;
って何だ?比較演算子として使ってるんぢゃなさそうだよね、この >=
配列data よりも大きな値をもつ dataEv の部分を配列 A に代入してるとか?

>「一家に一台・MATLAB」!
が無理ゲーなのだけは理解できました。キリッ
Commented by supernova1987a at 2020-01-06 16:55
> te kureさん
今思うと、昨年までの解析ごっこは相当に無茶な取り組みでした。今年の解析ごっこは昨年比100倍くらいのスピードで進行しそうな予感がします。
初めてパソコンを買った日、初めてワークステーションを与えられた日、などと同等の感動があります。
Commented by supernova1987a at 2020-01-06 17:50
> kem2017さん
行列変数がデフォという点にひかれてポチったものの、初日は難儀しました。
ちなみに、
C= A<Bと書くと、
AとBの同じ位置の要素を比較して、
BがAより大きければ1、そうでないなら0をCの同じ場所に書き込むという処理となります。
このCをBに掛け算させると、Bの要素のうちAよりも大きな物だけが残ります。
Commented by UTO at 2020-01-06 18:55 x
おぉ、ついに「まとらぼ」の呪文を覚えたのですね!

やー・・これ、自分もだいぶ前にちょっとだけ学ぶ機会があったのですが、結局、使わず終いになってしまい・・^^;
カメラ制御も上手くAPIがあればできちゃいますよね、、まとらぼ。

それにしても、素晴らしい解析ですね。
>一般的な「画像内の輝度揺らぎの標準偏差を測定する」方法では縮緬ノイズは評価できない
この辺りは、やっぱりユーザで実際に画像を扱ってるからこそですし、凄いなーと感心します。
あと、ダーク128枚あれば縮緬が減るとか、とても参考になります。
やっぱりCMOSは枚数勝負が一番かな~・・
Commented by 是空 at 2020-01-06 19:01 x
おいらは”MATLAB”よりも”MAD LOVE”の方が・・・

2020年の新しいオモチャを手に入れましたね。^^
しかしインタープリタなのに早いというのが不思議。
勝手な想像ですが、(←ここ重要)
解析計算って同じ計算をひたすら繰り返すというパターンが多いので、変数を配列にすることで計算のパターンを読み取れば(←インタープリタ)、あとはMATLAB内部だけでほとんど処理してしまうということなんでしょうね??
あと、マルチコア(スレッド)も効果的に使われそう。
Commented by えあ at 2020-01-07 00:09 x
すご!一気に解析が進みましたね。
面白い傾向も見えてきましたし、続きもたのしみです。
MATLABもすごいですけど、あっという間に使いこなしている あぷらなーとさん がすごすぎですよ。
Commented by kaz at 2020-01-07 04:12 x
学生時代使ってたが、貧乏なのでホームライセンスですら買えない、ってか買っても
Windowsが無いから使えないのだが「ラズパイ用」のMATLABはタダだそうで、これ、
ラズパイと言うかAstroberryでも動きませんかね?
100連装とは言わないから、4連装でLRGB一気に撮って4個のラズパイで処理してしまうとか。
Commented by rimpa at 2020-01-07 12:36 x
はじめてコメントします。
MATLAB互換のoctaveを使ってますが、fitsファイルを読み込ませたことはありませんでした。
確かに、行列の処理が簡単にできるので、複数画像のアライメントとかもできそうですね。
ありがとうございます。
Commented by M87JET at 2020-01-08 10:17 x
はじめまして。
M87JETと申します。
よくわかんないんですが、ものすごく未来を感じました。
ピクセル毎に特性を現場でマッピングして、
オーダーメード医療のように全素子キャリブレーションしてデータ処理して、
ハッブル宇宙望遠鏡超えの画像が得られるんじゃあないでしょうか?
茨城から西の方へ、お祈り。へへええ
Commented by noritos1047 at 2020-01-08 20:36 x
新年早々、何たる邪悪なブログであろうか。もうついていけん。
わしはもう短時間一枚撮り、ダークもフラットも省略じゃ。
と思ったのですが、短気を起こさずのんびり楽しむことを今年の抱負としました。
それにしてもこの消化不良、どうしてくれよう。パンシロンでも飲もうか。
Commented by supernova1987a at 2020-01-08 23:47
> UTOさん
お恥ずかしながら「まとらぼ」の存在を知りませんでした。まだ本領を発揮する機能は使っていませんが、自分がやりたかったことと相性が良くて驚愕です。

これまでモヤモヤしていた課題が一気に進み始めました。さて、何が出ますやら・・。(笑)
Commented by supernova1987a at 2020-01-09 08:27
> 是空さん
新年早々、良いオモチャを入手してテンションが上がりまくりです。おっしゃるとおり、通常のループ演算と行列変数を活かした演算とで非常に大きな差がありました。これまでできなかった諸々の課題が実現しそうです。
Commented by supernova1987a at 2020-01-09 08:31
> えあさん
本来、MATLABは高度なシミュレーションやディープラーニングなど多彩な機能を持つのですが、私にとってはFITSの読み書きと行列変数を用いた基礎計算が高速に実行できるだけで充分感激です。まだまだ少し触っただけですが、色んな面白い結果が出せそうです。
Commented by supernova1987a at 2020-01-09 08:36
> kazさん
ラズパイは触ったことがないので、動作するのか分かりませんが、ラズパイ4連装というアイデアは面白いですね。MATLAB自体はユーザーさんが多そうなので、すでに実現されている方がいらっしゃるかもしれません。
Commented by supernova1987a at 2020-01-09 08:40
> rimpaさん
ようこそ、拙ブログへ。octaveでFITSを扱うツールもあるみたいですね。おっしゃるとおり、行列変数の処理を上手く扱えば、アライメントのルーチンも書けそうです。私は使い始めたばかりなので、まだまだ道のりが遠いですが、いつか実現したいです。
Commented by supernova1987a at 2020-01-09 08:44
> M87JETさん
ようこそ、拙ブログへ。
アマチュア天文界ではピクセル個々の特性にこだわる酔狂な方は少なそうなので、逆に面白そうだと感じました。最終的に王道に帰着するかもしれませんが、ひょっとしたら何か新しい事が発掘されるかも?と期待しています。
Commented by supernova1987a at 2020-01-09 08:53
> noritos1047さん
ははは。新年早々、邪悪オーラ全開の記事になってしまいました。私は結果よりも過程を楽しむ方なので、正解にこだわらず色々と試してみたいと思います。あまりにも急に走ってしまったので書いてる方が消化不良に陥ってます。ご注意くださいませ。
名前
URL
削除用パスワード
by supernova1987a | 2020-01-06 08:07 | Comments(18)

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


by あぷらなーと