ほしぞloveログ

天体観測始めました。

2019年03月

ASI294MC Proのコンバージョンファクターを求める一環であぷらなーとさんからのコメントがあり、debayerせずにRGGBの画素ごとに分離してみたらという助言がありましたので試してみました。 

RGB分離したファイルからS-N^2分布をプロットするpythonコードを、fitsから直接CFA (Color Filter Array)のRGGBの4つ分に分離するようなものに改良して再度プロット。画像数は縦横2分の1、面積では4分の1になります。結果が以下のようになります。

20190304_02_Conversion_Factor_CFA



ところがグラフを見てもわかるように、全然メーカー値(コンバージョンファクター: 4程度、ユニティーゲイン:120程度)を再現できていません。要するに各画素の画像ファイルのノイズが大きい、すなわち輝度のばらつきが大きすぎるのです。 自分で書いたルーチンが何かおかしいのかと疑い、PixInsightを使ってマニュアルでも解析してみました。

PixInsightにはPreprocessingのところにSplitCFAというコマンドがあります。これを画像に適用してやると、画素ごとに分離された画像が自動的に出来上がります。この画像を使って、これもPixInsight上のImageInspection -> Statisticsで統計的な値をみてやりましたが、pythonで解析したのと同様の傾向でした。また、CosmeticCorrectionでのホットピクセルとクールピクセルの除去は結果にほとんど影響を与えませんでした。これは中心付近だけをみてもほとんどホットピクセルとクールピクセルが入っていないからかと思われます。

わかりやすいところで、輝度10000付近になるところを例として、結果を示しておきます。

Capture_00_16_30_00010_00_20_06_cc_CFA0->Preview01
CFA0 CFA1 CFA2 CFA3
count (px) 15544 12540 15120 14868
mean 10205.3 9805.5 9803.4 9209.9
median 10205.5 9806.6 9803.6 9209.6
stdDev 65.7 67.0 65.3 61.1
avgDev 52.6 53.2 51.9 48.9
MAD 44.0 44.0 44.0 42.0
minimum        9912.5 9528.6 9543.6 8977.6
maximum        10444.5 10049.5 10084.5 9475.6


Capture_00_16_30_00010_00_20_06_cc_RGB_VNG->Preview01
R G B
count (px) 14950 14950 14950
mean 10229.3 9798.6 9206.1
median 10229.5 9798.8 9206.9
stdDev 59.566 56.9 57.6
avgDev 47.138 44.4 45.6
MAD 39.040 36.2 38.0
minimum 10003.5 9551.6 8941.1
maximum 10484.4 10049.5 9445.2


上がCFAで分離した場合、下がdebayerした場合になります。赤いところを見比べればいいのですが、全く同じ画像から作り出したにも関わらず、明らかにCFA分離の方が平均偏差(avgDev)が大きく出てしまっています。あまり大きな差に見えないかもしれませんが、分布図を作るときにこの値を2乗してプロットするので、ざっくり(52/46)^2で1.28、約3割の増加です。1次グラフの傾きも3割くらい変わってきて、コンバージョンファクターも、ユニティーゲインも3割くらい(ZWOカメラのゲインにして40程度)変化があります。

どうもPixInsightでdebayerをするときにばらつきが減るような効果があるようなのですが、なぜそうなるのかはまだよく理解できていません。ただ、SharpCapでもただちにRGBにdebayerして測定していると思われるので、とりあえずはdebayerを正しいと思って進めることにします。

まとめですが、今回の過程でまだ2つのしっくりこない部分が残りました。
  • ばらつき具合を評価するのに標準偏差よりも平均偏差を使ったこと。
  • debayerとCFA分離でノイズに明らかに差が出る。debayerの方がノイズが少ない。
ということが判明したのですが、特に後者のわかりやすい説明があれば知りたいです。

debayerの時に周りの(同じカラーの)ピクセルの値を一部情報として輝度を決めるはずなので、確かに少しばらつきがなくなるというのも一理あるのですが、はたしてそれを性能評価として使っていいものなのか?科学的には間違っている気がしてなりません。

ではなぜメーカー値もSharpCapも同じ値を出せるのか?やはりまだ何か自分が根本的に間違っているのでしょうか?


昨日の記事の続きです。やっとなんとか形になりました。

ノイズを標準偏差で評価するか、平均偏差で評価するか迷っていたのですが、Twitteでガウス分布から外れているのなら飛んだ値が多いはずなので、(飛んだ値に影響されにくい)平均偏差の方がいいのではという意見をもらいました。なるほど、考えてみればその通りで、標準偏差と平均偏差にすでに無視できないような有意な差があるということは、いいかえてみればガウス分布から外れた値も多いということが言えるのではと思います。

Fits画像のhistogram


というわけで実際にヒストグラムで分布を見てみました。まず、debayerなど何の処理もしていないRAW画像です。 

histgram_raw

見ての通り、ガウス分布からかなり外れていることがわかります。これはRGBでそれぞれ反応が違うために山がいくつもできるのかと思われます。

次に、同じ画像をdebayerしてRed、Green、Blueに分けたヒストグラムものを示します。

まずはRed:
histgram_R
次にGreen:
histgram_G

最後にBlue:
histgram_B


不思議なのは、RGBに分離しても山がいくつも見えることです。debayerの際に周りのピクセルの状況も読み込んでいるからなのか、もしくは画面の中で場所によって明るさに違いがあって、それが山になっているのかもしれません。また、RGBを合わせてもRAWの山の形になりそうもないことも不思議です。一瞬違う画像を処理したかと思ってしまったのですが、きちんと確認しても同じRAW画像から分離したものです。debayerもそんなに単純でないようです。

最後に、その中の50x50ピクセルを取り出してきた場合のヒストグラムです。
histgram_50x50
山がいくつもあるようなことはなくなり、大まかな形としてはガウス分布にだいぶん近づきます。それでもサンプル数が少ないことによるばらつきがあるのも確かなので、ここでは平均偏差でいくのが良いと考えることにします。



Conversion Factor

さて、実際にコンバージョンファクターを求めてみました。サンプル数を多くするために画像中心付近の100x100ピクセルを選んで解析しています。

結局今回はPythonで平均偏差を求めるルーチンをを自前で書いて、各ピクセルごとに計算しています。書き忘れてましたが上のヒストグラムも全部Pythonで書いています。やっとPythonでの画像解析に慣れてきました。結果ですが、以下のようになります。

20190302_01_Conversion_Factor

ついにここまでくることができました。結果はグラフの中にも数値で書いてありますが、コンバージョンファクターとして4.12、そこから計算できるUnity gainが200 x log10(4.12) = 123となり、メーカー値の117とわずか0.6dB、1.07倍の誤差くらいの範囲で求めることができました。

検証


もう少し検証してみます。

IMG_3262

上のようなSharpCapでの自動測定の結果のグラフと比べると、自動測定の測定値を伸ばしていくと0点近くに行きますが、自分で測定したものは0点に向かわずに、y切片で-352くらいのところにあたります。本当にきちんと測定しようとするならバイアスノイズをのぞいたり、フラット補正をすべきなのですが、今回は省いています。それでもSharpCapもそれ専用の測定はしていないように見えるので、うまくy切片が0になるような補正をかけているものと考えられます。

もう一点、自動測定の場合、測定点がいくつか重なっているように見えます。おそらくこれはRGBと分解した3点が重なっていると推測されるのですが、それにしても横軸(ADU)が一致しすぎています。普通に測定すると、自分で測定した時みたいにRGBで光源も違えばセンサーのフィルター特性も違うはずなので、ずれるはずです。これもSharpCapの自動測定では何らかの補正をしているものと思われます。


まとめ

結局、上の結果を得るまでに2週間くらいかかりました。色々苦労しましたが得たものも多く、まずPythonでの画像解析の環境がだいぶ揃いました。既存ライブラリに頼らない、ピクセルごとに解析する手法もある程度得ることもできました。統計的にどのようにアプローチすればいいのかも少し学ぶことができました。

次はEOS 6Dのユニティーゲインを求めることでしょうか。
あー、ホントはCP+行きたかったです。

今回はCMOSカメラ、ZWOのASI294MC Proの性能評価の一環で、全ての測定の元になるADUからeへの変換のコンバージョンファクターの測定についてです。結論から言いますと、SharpCapの自動測定機能での結果と、SharpCapでマニュアルで一枚一枚撮影しその画像を自分で解析した結果がどうしても合いません

この記事は多分ほとんどの人にはめんどくさい話で、よほどでない限り興味がないことと思いますし、しかもうまく結果が出なかったものなので、公表するかどうかも迷っていたのですが、それでも自分のメモがわりに書いておこうと思います。ご容赦ください。


動機

もともとダークノイズを評価する過程の一環で進めているのですが、今回の測定の動機は2つあります。
  1. ダークノイズの測定は多岐に渡るので、まずは解析環境をpythonで整えようとするのにちょうどいい練習になる
  2. 天体用CMOSカメラだけでなく、一眼レフカメラの性能評価もできないかと思ったから
です。特にEOS 6Dのユニティーゲイン(ADUとeの比が1になるゲイン)、引いてはコンバージョンファクターの測定まで自前でできたらなと目論んでいたのですが、今のところ見事に失敗。

最近ブログをなかなか更新できなかったのは、天気が悪くて星が見えないとか、仕事が忙しいとかもありますが、この解析が全然うまくいかなくてずっと悩んでいたというのが、一番大きな理由です。


測定方法

各画像の撮影はSharpCapで撮影します。共通の設定は
  • iPadのColor Screenというソフトを光源とした
  • RAW16
  • Gain = 0
  • Brightness = 8
  • White Bal(R) = 50
  • White Bal(B) = 50
  • 温度15度程度(コントロールなし)
となります。この状態で露光時間を変更して10枚程度の画像を撮影します。上記設定や露光時間はSharpCapでのセンサー性能を測る時のパラメーターを参考にしています。というか、最初適当に設定していたのですが、結果が全然合わないので、最後はコンバージョンファクターを測る時の状況に限りなく合わせるようにしました。

ゴールとしては下の写真(SharpCapのセンサー性能測定機能で自動測定した場合)の

IMG_3260

右のようなグラフが得られればOKです。横軸(各ピクセルの明るさ)が10000程度の時に縦軸(ノイズの2乗)が2500程度です。グラフの傾きは0.25程度、その傾きの逆数が今回求めたいコンバージョンファクターになり、普通に測定すると1/0.25=4程度になるはずです。この値はZWOが示している値ともほぼ一致しています。

コンバージョンファクターはちょっと理解が大変かもしれませんが、関係式と意味についてはこのページの1のところに、式の証明についてはこのページの一番最後のおまけのところに書いてあります。


測定結果

ところが自分で測定してみた結果は散々なものです。普通に画像を撮影してそのまま何も考えずに解析すると、そもそもDebayerもされていなかったりするので、ノイズが大きく出すぎてしまいます。結果を見せるのもあほらしいのですが、

mymag_all


のようになり、SharpCapの結果にカスリもしないくらいノイズが大きく出てしまっています。傾きが30くらい、コンバージョンファクターは0.03とかで、メーカー値の100分の1以下です。ここから苦難のノイズハンティングの道が始まりました。結局やったことをまとめると
  1. SharpCapでfitsファイルを撮影
  2. PixInsightでCosmeticCorrectionでホット、クールピクセルを除去 (飛び抜けて明るいピクセルなどあるとばらつきが大きく出てしまう)
  3. PixInsightでDebayerをしてカラー化 (RGBでゲインが多少違うため、debayerせずに標準偏差を取るとばらつきが大きく出てしまう)
  4. PixInsightでR、G、B画像に分離し、一枚一枚を個別に保存する (今回解析に使ったirafは天文研究に使われるソフトで、モノクロがほぼ前提なので、カラー画像を解析できない)
  5. 中心近くの50x50ピクセルのみを選択して解析 (画像全体だと周辺減光などの影響で、ばらつきが大きく出てしまう)

4番まで進めた時のグラフが

mymag _rgb_all


のようになりますが、まだ傾きが10程度、コンバージョンファクターにして0.1程度しかありません。

さらに5番目の中心部分のみを解析するようにして、やっと下のグラフくらいにまでなりました。

mymag_rgb_cut


それでも結局傾きが0.35程度、コンバージョンファクターが1/0.35で3程度になり、どうしてもまだ4近くにまでなりません。

考察と今後

なぜこの差が縮まらないか、もう少し検証します。まず、横軸(各ピクセルの明るさ)が10000の時に縦軸(ノイズの2乗)が2500くらいになるためには、ノイズはそのルートの50程度でなければなりません。ではノイズと言っているものが何かと言うと、画像から測定したピクセルの明るさのばらつきということなので、普通は標準偏差(standard deviation)をとればいいと思われます。この標準偏差を求めるのに今回は天文研究でよく使われているirafを使いました。ところが、明るさ10000程度の50x50ピクセルの明るさのばらつきの標準偏差をirafで測定すると60程度になってしまいます。グラフの横軸でいうと2乗なので3600程度になってしまうわけです。

ここでirafを疑いました。何か間違った結果が出ているのではと思ったのです。そこでPixInsightの統計ツールで測定したのですが、標準偏差はやはり60程度とでます。それどころか、SharpCapでも画像の選択したあるエリアの各色の標準偏差をリアルタイムで測定できるのですが、それもやはり60程度なのです。

SharpCapで測定しても60とでるならば、SharpCapの自動センサー性能測定の測定はどうやってやっているのでしょうか?何か特別なことをやって50と出しているのか、それともまだ私が何か勘違いをしているのか

PixInsightの統計ツールで少しヒントになるようなものを見つけました。標準偏差ではなくてオプションでAverage absolute deviationという値を出すことができるのですが、この値がちょうど50程度になるのです。

IMG_6410
標準偏差(stdDev)が60ちょい、Average absolute deviation(avgDev)が
50切るくらいになっているのがわかると思います。

Average absolute deviationのは一般的にはMean absolute average (around the mean)というらしくて日本語では単純に平均偏差というらしいです。標準偏差が各値(Xi)から平均値(M)を引いたものを2乗したものの総和を総数Nで割ったもの

1Ni=0N(XiM)2

に対して平均偏差は各値から平均値を引いたものの絶対値の総和を総数で割ったもの

1Ni=0N(XiM)

となります。他にもMean absolute average (around the median)というのもありますが、こちらは平均値を引く代わりに中心値を引きます。

標準偏差が2乗和のルートになるので、ばらつきがより効いてくることになり一般的に

標準偏差 > 平均偏差

となるそうで、確かに標準偏差より小さい値になっていて納得です。

さて、平均偏差を使えば、メーカー値もしくはSharpCapで測定した値に近い結果が出るはずなのですが、そもそも平均偏差を使っていいものなのか?やはり普通に考えると標準偏差を使った方が、あとの統計的な評価が簡単になりそうで、素直な気がします。

さらに、irafなどの一般的な解析ツールでは平均偏差を出すことができるものが少ないので、グラフまで出せるくらいにきちんと解析するのなら自前で統計処理の部分のコードを書く必要があります。

そんなこんなで、今pythondで書いているのですが、果たしてこの方向が正しいのかどうか?
まだ色々迷っています。


このページのトップヘ