ほしぞloveログ

天体観測始めました。

タグ:コンバージョンファクター

2022年の反省でも述べましたが、最近自分で考えることがあまりできてなかったので、今回の記事は久しぶりに計算です。内容は、撮影した画像にはどんなノイズが入っていて、それぞれどれくらいの割合になっているかを見積ってみたという話です。


動機

まずは動機です。1月に開田高原でM81を撮影した時に、M81本体に加えて背景のIFNが見えてきました。下は300秒露光を28枚撮影したL画像を強度にオートストレッチしています。
masterLight_BIN-2_4144x2822_EXPOSURE-300.00s_FILTER-L_mono

一方、2022年の5月にも自宅で同じM81を撮影しています。背景が埃とスカイノイズで淡いところがまってく出なかったので記事にしていないのですが、下は露光時間600秒を22枚撮影したL画像で、強度にオートストレッチして、かつできるだけ見やすいようにABEとDBEをかけてカブリを排除しています。開田高原の時よりもトータル露光時間は1.6倍ほど長く、被りを除去しても、淡い背景については上の画像には全く及びません。
masterLight_600_00s_FILTER_L_mono_integration_ABE_DBE3
明らかにスカイノイズの影響が大きいのですが、これを定量的に評価してみたくなったというものです。

もう一つの動機ですが、このブログでも電視観望について多くの記事を書いています。電視観望はリアルタイム性を求めることもあるので、露光時間を短くしてゲインを高く設定することが多いです。この設定が果たして理に適っているのかどうか、これも定量的に議論してみたいとずっと思っていました。


目的

これからする一連の議論の目的ですが、
  1. 画像に存在するどのノイズが支配的かを知ること。
  2. 信号がノイズと比較して、どの程度の割合で効くのかを示す。
  3. 電視観望で高いゲインが有利なことを示す。
を考えてたいと思っています。今回の記事はまずは1番です。こちらはスカイノイズをどう評価するかが鍵なのかと思っています。

2番は意外に難しそうです。信号である天体は恒星ならまだしも、広がっている天体の明るさの評価ってなかなか大変です。これは後回しにするかもしれません。

3番は長年の電視観望がなぜ短時間で天体を炙り出せるのかという疑問を、定量的に表す事ができたらと考えています。こちらは大体目処がついてきたので、近いうちに記事にするかと思います。


基本的なノイズ

天体画像撮影におけるノイズは5年くらい前にここで議論しています。SN比の式だけ抜き出してくると

S/N=ntSsigAnσ2+AntSsky+AntSdark+ntSsignS/N=ntSsigAnσ2+AntSsky+AntSdark+ntSsign
と書くことができます。ここで、
  • AA [pix] : 開口面積
  • nn : フレーム枚数
  • σσ [e-/pix] : 読み出しノイズ
  • tt [s] : 1フレームの積分(露光)時間
  • SskySsky  [e-/s/pix] : スカイバックグラウンド
  • SdarkSdark  [e-/s/pix] : 暗電流
  • SsigSsig  [e-/s] : 天体からの信号
となり、S/Nとしては何枚撮影したかのルートに比例する事がわかります。今回のノイズ源としては読み出しノイズ、スカイノイズ、ダークノイズを考えます。ショットノイズは天体などの信号があった場合に加わるノイズですが、今回は天体部分は無視し背景光のみを考えるため、天体からのショットノイズは無視することとします。

重要なことは、読み出しノイズは露光時間に関係なく出てくるために分母のルートの中にtがかかっていないことです。そのため、他のノイズは1枚あたりの露光時間を伸ばすとS/Nが上がりますが、読み出しノイズだけは1枚あたり露光時間を増やしてもS/Nが上がらないということを意識しておいた方がいいでしょう。その一方、撮影枚数を稼ぐことは全てのノイズに対してS/N改善につながり、読み出しノイズも含めて撮影枚数のルートでS/Nがよくなります。繰り返しになりますが、一枚あたりの露光時間を伸ばすのでは読み出しノイズだけは改善しないので、他のノイズに比べて効率が悪いということです。

分子にあたる天体信号の評価は意外に大変だったりするので、今回は分母のノイズのみを考えることにします。


パラメータ

上の式を元に、ここで考えるべきパラメータを固定しやすいもの順に書いておきます。
  1. カメラのゲイン(gain)
  2. 温度 (temperature)
  3. 1枚あたりの露光時間 (time)
  4. 空の明るさ
の4つで実際に撮影した画像の各種ノイズを推定できるはずです。少し詳しく書いておくと、
  • 1は撮影時のカメラの設定で決める値です。読み出しノイズ (read noise)を決定するパラメータです。
  • 2は撮影時のセンサーの温度で、冷却している場合はその冷却温度になります。単位は[℃]となります。この温度と次の露光時間からダークノイズを決定します。
  • 3は撮影時の画像1枚あたりの露光時間で、単位は秒 [s]。2の温度と共にダークノイズ (dark noise)を決定します。ダークノイズの単位は電荷/秒 [e/s]。これは[e/s/pixel]と書かれることもありますが、ここでは1ピクセルあたりのダークノイズを考えることにします。なお、ホットピクセルはダークノイズとは別と考え、今回は考慮しないこととします。ホットピクセルは一般的にはダーク補正で除去できる。
  • 4は撮影場所に大きく依存します。今回は実際に撮影した画像から明るさを測定することにします。単位は [ADU]で、ここからスカイノイズを推測します。ちなみに、3の露光時間も空の明るさに関係していて、長く撮影すれば明るく写り、スカイノイズも増えることになります。

実際のパラーメータですが、今回の記事ではまずはいつも使っている典型的な例で試しに見積もってみます。私はカメラは主にASI294MM Proを使っていて、最近の撮影ではほとんど以下の値を使っています。
  1. 120
  2. -10℃
  3. 300秒
これらの値と実際の画像から背景光を見積もり、各種ノイズを求めることにします。


読み出しノイズ

読み出しノイズはカメラのゲインから決まります。ZWOのASI294MM Proのページを見てみると、


真ん中の少し手前あたりにグラフがいくつか示してあります。各グラフの詳しい説明は、必要ならば以前の記事



をお読みください。上の記事でもあるように、カメラの各種特性はSharpCapを使うと自分で測定する事ができます。実測値はメーカーの値とかなり近いものが得られます。

グラフから読み出しノイズを読み取ります。gain 120の場合おおよそ

1.8 [e rms]

ということがわかります。rmsはroot mean sqareの意味で日本語だと実効値、時系列の波形の面積を積分したようなものになります。例えば片側振幅1のサイン波なら1/√2で約0.7になります。他のノイズも実効値と考えればいいはずなので、ここではrmsはあえて書かなくて

1.8 [e]

としていいでしょう。

読み出しノイズは、実際の測定では真っ暗にして最短露光時間で撮影したバイアスフレームのノイズの実測値に一致します。以前測定した結果があるので、興味のある方はこちらをご覧ください。



ダークノイズ

ダークノイズの元になる暗電流に関しては温度と1枚あたりの露光時間が決まってしまえば一意に決まってしまい、これもZWOのASI294MM Proのページから読み取ることができます。

私はこの値を実測したことはないのですが、そーなのかーさんなどが実測していて、メーカー値とほぼ同じような結果を得ています。

グラフから温度-10℃のところの値を読み取ると暗電流は

0.007 [e/s/pixel]

となるので、1枚あたりの露光時間300秒をかけると

2.1 [e/pix]

となります。/pixはピクセルあたりという意味なので、ここは略してしまって

2.1[e]

としてしまえばいいでしょう。単位[e]で考えた時に、暗電流のルートがダークノイズになるので、

sqrt(2.1)=1.5[e]

がダークノイズとなります。

(追記: 2023//4/19) ちょっと脱線ですが、だいこもんさんのブログを読んでいると、2019年末の記事に「dark currentはgain(dB)に依存しないのか?」という疑問が書かれています。答えだけ言うと、[e]で見ている限り横軸のゲインに依存しないというのが正しいです。もしダークカレント、もしくはダークノイズを[ADU]で見ると、元々あったノイズがアンプで増幅されると言うことなので、単純に考えて横軸のゲイン倍されたものになります。実際の画面でも横軸ゲインが高いほど多くのダークノイズが見られるでしょう。でも、コンバージョンファクターをかけて[ADU]から[e]に変換する際に、コンバージョンファクターが横軸のゲイン分の1に比例しているので、積はゲインによらずに一定になるということです。

ちなみに、ホットピクセルはダーク補正で取り除かれるべきもので、ここで議論しているダークノイズとは別物ということを補足しておきたいと思います。

(さらに追記:2023/10/15)
ダークファイルを撮影したので、実際のダークノイズを測定してみました。画像からのノイズの読み取り値は2.6 [ADU]でした。コンバージョンファクターで単位を[e]にすると、2.6x0.9 = 2.3 [e]。ここから読み出しノイズを引いてやると、2乗の差のルートであることにちゅいして、sqrt(2.3^2-1.8^2) = 1.5 [e] と見積り値に一致します。このように、見積もりと実測が、かなりの精度で一致することがわかります。


スカイノイズ

ここが今回の記事の最大のポイントです。読み出しノイズとダークノイズだけならグラフを使えばすぐに求まりますが、恐らくスカイノイズの評価が難しかったので、これまでほとんど実画像に対してノイズ源の評価がなされてこなかったのではないでしょうか?ここでは実画像の背景部の明るさからスカイノイズを推測することにします。

まず、明るさとノイズの関係ですが、ここではコンバージョンファクターを使います。コンバージョンファクターはカメラのデータとして載っています。例えばASI294MM Proでは先ほどのZWOのページにいくと縦軸「gain(e/ADC)」というところにあたります。コンバージョンファクターの詳しい説明は先ほど紹介した過去記事を読んでみてください。コンバージョンファクターは他に「ゲイン」とか「システムゲイン」などとも呼ばれたりするようです。名前はまあどうでもいいのですが、ここではこの値がどのようにして求められるかを理解すると、なぜスカイノイズに応用できるか理解してもらえるのかと思います。

コンバージョンファクターの求め方の証明は過去記事の最後に書いてあるので、そこに譲るとして、重要なことは、画像の明るさSとノイズNは次のような関係にあり、
(N[ADU])2=S[ADU]fc[e/ADU](N[ADU])2=S[ADU]fc[e/ADU]明るさとノイズの二つを結ぶのがコンバージョンファクターfcとなるということです。逆にいうと、このコンバージョンファクターを知っていれば、明るさからノイズの評価が共にADU単位で可能になります。

もっと具体的にいうと、コンバージョンファクターがわかっていると、スカイノイズが支配していると思われる背景部分の明るさを画像から読み取ることで、スカイノイズを直接計算できるということです。これは結構凄いことだと思いませんか?


実際のスカイノイズの見積もり

それでは実画像からスカイノイズを見積もってみましょう。最初に示した2枚のM81の画像のうち、上の開田高原で撮影した画像の元の1枚撮りのRAW画像を使ってみます。

上の画像は既にオートストレッチしてあるので明るく見えますが、ストレッチ前の実際のRAW画像はもっと暗く見えます。明るさ測定はPixInsight (PI)を使います。PIには画像を解析するツールが豊富で、今回はImageInspectionのうち「Statistics」を使います。まず画像の中の暗く背景と思われる部分(恒星なども含まれないように)をPreviewで選び、StatisticsでそのPreviewを選択します。その際注意することは、カメラのADCのbit深度に応じてStatisticsでの単位を正しく選ぶことです。今回使ったカメラはASI294MM Proなので14bitを選択します。輝度の値は「mean」を見ればいいでしょう。ここでは背景と思われる場所の明るさは約920 [ADU]と読み取ることができました。ついでに同じStatisticsツールのノイズの値avgDevをみると17.8 [ADU]と読み取ることが出来ました。

もっと簡単には、画像上でマウスを左クリックするとさまざまな値が出てきますので、その中のKの値を読み取っても構いません。ここでも同様に、単位を「Integer Renge」で手持ちのカメラに合わせて「14bit」などにすることに注意です。

いずれのツールを使っても、背景と思われる場所の明るさは約920[ADU]と読み取ることができました。

前節の式から、輝度をコンバージョンファクターで割ったものがノイズの2乗になることがわかります。gain120の時のコンバージョンファクターはグラフから読み取ると0.90程度となります。

これらのことから、背景に相当する部分のノイズは以下のように計算でき、

sqrt(920/0.90) = 32.0 [ADU] -> 28.8 [e]

となります。どうやら他の読み出しノイズやダークノイズより10倍程大きいことになります。あれ?でもこれだとちょっと大きすぎる気がします。しかも先ほどのStatisticsツールでの画面を直接見たノイズ17.8[ADU]をコンバージョンファクターを使って変換した

17.8 [ADU] x 0.90 [e/ADU] = 16.0 [e]

よりはるかに大きいです。これだと矛盾してしまうので何か見落としているようです。

計算をじっくり見直してみると、どうやら測定した輝度は「背景光」とオフセットの和になっているのに気づきました。撮影時のオフセットとして40を加えてありますが、この値に16をかけた640がADUで数えたオフセット量として加わっているはずです。実際のマスターバイアス画像を測定してみると、輝度として平均で約640 [ADU]のオフセットがあることがわかったので、これは撮影時に設定したものとぴったりです。この値をを920 [ADU]から引いて、280 [ADU]を背景光の貢献分とします。その背景光からのスカイノイズ成分は

sqrt(280/0.90) = 17.6 [ADU]


これを[ADU]から[e]に変換するためにさらにコンバージョンファクター[e/ADU]をかけて

17.6 [ADU] x 0.90 [e/ADU] = 15.9 [e]

となります。これだと画面からの実測値16.0[e]と少なくとも矛盾はしませんが、既に実測のトータルノイズにかなり近い値が出てしまっています。果たしてこれでいいのでしょうか?


各ノイズの貢献度合い

次にこれまで結果から、各ノイズの貢献度を見ていき、トータルノイズがどれくらいになるのかを計算します。

画像1枚あたりの背景に当たる部分の各ノイズの貢献度は多い順に
  1. スカイノイズ: 15.9[e]
  2. 読み出しノイズ: 1.8 [e]
  3. ダークノイズ: 1.5[e]
となります。Statisticsツールで測定した実際の1枚画像の背景のノイズの実測値は14.9[e]程度だったので、上の3つのノイズがランダムだとして2乗和のルートをとると、

sqrt(15.9^2+1.8^2+1.5^2) = 16.0 [e]

となり、実測の16.0[e]になる事がわかります。スカイノイズに比べてトータルノイズがほとんど増えていないので、読み出しノイズとダークノイズがほとんど影響していないじゃないかと思う方もいるかもしれません。ですが、互いに無相関なノイズは統計上は2乗和のルートになるので本当にこの程度しか貢献せず、実際にはスカイノイズが支配的な成分となっています。

ここまでの結果で、今回のスカイノイズ成分の推定は、定量的にも実画像からの測定結果とほぼ矛盾ないことがわかります。この評価は結構衝撃的で、暗いと思われた開田高原でさえスカイノイズが圧倒的だという事がわかります。


富山の明るい空

ちなみに、最初の2枚の画像のうち、下のものは自宅で撮影したM81です。富山の明るい北の空ということもあり、そのRAW画像のうちの最も暗い1枚(2022/5/31/22:48)でさえ、背景の明るさの読み取り値はなんと3200[ADU]程度にもなります。バイアス640を引くと2560[ADU]で、開田高原の背景光の値240に比べ約10倍明るく、ノイズは

sqrt(2560/0.90) = 53.3 [ADU] -> 48.0 [e]

となります。実際の画像からPIのStatisticsで読み取ったトータルノイズの値が53.4[ADU]->48.1[e]でした。

スカイノイズ48.0[e]に、読み出しノイズ1.8 [e] 、ダークノイズ1.5[e]を2乗和のルートで考えるとトータルノイズは

sqrt(45.0^2+1.8^2+1.5^2) = 48.1 [e]

となり、明るい画像でも背景光の輝度から推測したノイズをもとに計算したトータルノイズ値と、画面から実測したノイズ値が見事に一致します。

開田高原の暗い画像でさえスカイノイズが支配的でしたが、富山の明るい空ではスカイノイズが読み出しノイズ、ダークノイズに比べて遥かに支配的な状況になります。淡いところが出なかったことも納得の結果です。

うーん、でもこんなスカイノイズが大きい状況なら
  • もっと露光時間を短くして読み出しノイズを大きくしても影響ないはずだし、
  • 冷却温度ももっと上げてダークノイズが大きくなっても全然いいのでは
と思ってしまいます。でもそこはちょっと冷静になって考えます。早急な結論は禁物です。次回の記事では、そこらへんのパラメータをいじったらどうなるかなどに言及してみたいと思います。


まとめ

背景の明るさとコンバージョンファクターから、スカイノイズを見積もるという手法を考えてみました。実際の撮影画像のノイズ成分をなかなか個別に考えられなかったのは、スカイノイズの評価が難しいからだったのではないかと思います。

今回の手法は背景の明るさが支配的と思われる部分がある限り(天体写真のほとんどのケースは該当すると思われます)スカイノイズを見積もる事ができ、状況が一変するのではと思います。また手法も輝度を測るだけと簡単なので、応用範囲も広いかと思われます。

今回の手法を適用してみた結果、実際に遠征した開田高原のそこそこ暗い空で撮影した画像でさえもスカイノイズが支配的になる事がわかりました。もっと暗い状況を求めるべきなのか?それとも露光時間を短くしたり、温度を上げてしまってもいいのか?今後議論していきたいと思います。

とりあえず思いつくアイデアをばっとまとめてみましたが、もしかしたら何か大きな勘違いなどあるかもしれません。何か気づいたことがありましたらコメントなどいただけるとありがたいです。


この記事の続き



になります。露光時間と温度について、実際の撮影においてどの程度にすればいいか評価しています。









前回の6Dのユニティーゲインの記事ですが、難しいという話を聞いたので、できるだけわかりやすく解説してみようと思います。




コンバージョンファクター


IMG_2032

何はともあれ、まず重要なのはコンバージョンファクター (conversion factor) とか、gain (ややこしいのですがISOに相当するgainとは全然の別の意味です。以下区別するために、コンバージョンファクターの意味ではgainやゲインという言葉は使わず、ISOの意味でのみgainもしくはゲインと書きます。)とか呼ばれている、センサーで発生する電子の数と出力信号のカウント数を変換する係数で、単位は[e/ADU]になります。eは電子1個の単位、ADUはADC (Analog to Digital Converter) でカウントする1カウントの単位です。発生する電子の数は後で書きますが、検出される光子の数と比例するので、ここではとりあえず光子の数と言ってしまうことにします。

このページの結果によるとEOS 6Dの場合、例えば、ISO100のとき、コンバージョンファクターは5.6413 [e/ADU]で、光子が6個近くセンサーの1素子に入ってやっとADCのカウントが1増えます。6Dの場合14bit = 16384なので、5.6413 x 16384 = 92427個の光子が1素子に入ると、その素子のカウントは一杯になり「飽和」状態になります。このブログでは「サチる」とか、「サチった」とかいう表現をしています。これは「飽和」の英語Saturationから来ています。

例えば、ISO400のとき、コンバージョンファクターは1.4178 [e/ADU]となり、光子が1.5個くらい入るとADCのカウントが1増えます。

この表から考えると、ISO575くらいが実現できるなら、コンバージョンファクターは1 [e/ADU]となり、光子が1個入るとADCのカウントが1増えます。このときのゲインをその名の通り、ユニティーゲイン(unity gain)と呼びます。ユニティーは1という意味ですね。ISOなので、ユニティーISOとか読んでもいいでしょう。呼び方はまあどうでも良くて、重要なのはセンサーで検出される光子1個がADCを1カウント増やすという、1対1の関係です。

CMOSセンサーの解析はこのコンバージョンファクターの値を求めるところから全てが始まります。


コンバージョンファクターの原理

ではコンバージョンファクターを求めるためにはどうしたらいいのでしょうか?まずは原理式を理解してみましょう。CMOSセンサーをある出力ゲインに固定して測定した信号\(S\mathrm{[ADU]}\)とそのときのノイズ\(N\mathrm{[ADU]}\)には以下の関係があります。

\[(N\mathrm{[ADU]})^2=\frac{S\mathrm{[ADU]}}{f_{\mathrm{c}}\mathrm{[e-/ADU]}}+\frac{(N_{\mathrm{read}}\mathrm{[e-]})^2}{(f_{\mathrm{c}}\mathrm{[e-/ADU}])^2}\]

このとき\(N_{\mathrm{read}}\mathrm{[e-]}\)は読み出しノイズ、\(f_{\mathrm{c}}\mathrm{[e-/ADU]}\)がコンバージョンファクターです。

右辺2項目の読み出しノイズは十分小さい仮定として、簡単に

\[(N\mathrm{[ADU]})^2=\frac{S\mathrm{[ADU]}}{f_{\mathrm{c}}\mathrm{[e-/ADU]}}\]
を考えます。

画像を撮影して、その1ピクセルの明るさ\(S\mathrm{[ADU]}\)を測り、その1ピクセルの明るさがどれくらいバラけているかを多数プロットしてやればいいのです。

この式自身の証明はこのページの最後のおまけの部分を見てください。ちょっととっつきにくいと思うかもしれませんが、ショットノイズの関係式だけから数学的に綺麗に出てくるので、話としては至極単純です。逆にこの関係式があるので、多数の点数をとってくれば統計的にコンバージョンファクターが確定するというわけです。多数の点をとってくるのは、画像ファイルが多数の点でできていることを考えると十分可能で、アイデア次第で多くのサンプルを取り出すことができるわけです。この関係式をものすごくうまく利用してますよね。

多くのサンプルを取り出すのはいろいろな方法があります。
  1. 一画面内に暗いところから明るいところまで写っている、同じ画角の景色などを多数枚撮影し、ある位置のピクセルに注目し、その平均値とバラけ具合を多数枚にわたって見る。多数のピクセルに対して同じことをする。
  2. フラットに近い画像を露光時間を変えて何枚か撮り、一枚のある100x100くらいのエリアの平均値とそのバラけ具合を見る。露光時間ごとにプロットする。
  3. 星などの適当な画像を2枚撮影し、その2枚の画像の差を取り(信号を消して、ノイズの2乗和のルートをとるということ)、その中の明るさが一定のエリアの平均値とバラけ具合を見る。明るさの違う領域で、同じことをしてプロットする。
などがあります。私が一番最初に学んだ方法は1.でした。SharpCapは2.の方法をとっています。最初に紹介したページ(大元はこのページです)やPixInsightでは3.を使っています。3.が撮影枚数が2枚と少なく、一番簡単でしょうか。工夫次第でまだ測定方法はいろいろ考えることができると思います。

PixInsightでの測定方法はNiwaさんが秀逸なタイトルをつけて詳しく解説してくれています。




実際の測定例

実際に信号とそのバラつきをプロットしたグラフを見てみましょう。まずは2.の例のSharpCapで測定した場合です。6Dの計測の時にグラフを写真に取り損なったので、ASI294MCの測定の時の写真を示します。

IMG_3262

一番右が、横軸明るさS、縦軸ノイズNの2乗でプロットしたものになります。測定点を結ぶと一本の線になり、その傾きの逆数がコンバージョンファクターになります。この場合、横軸が5000くらいの時に縦が1300くらいなので、傾きは0.26、その逆数は1/0.26=3.85[e/ADU]くらいになります。すなわち、光子3.85個入って、やっとADCのカウントが1進むということです。

次の例は自分で画像を撮影して測定した時の結果です。SharpCapがやる過程をマニュアルでやったような形になります。すなわち、同じゲインで露光時間を変えて何枚か撮影し、あるエリアの明るさとバラけ具合を測定すると言うものです。画像解析はMaltabを使ってやりました。Matlabを使ったのは、画像読み込みが楽なのと、平均や分散などの統計解析が揃っているからです。別に一枚一枚Photoshopとかで解析しても原理的にはできるはずです。センサーはASI290MMでモノクロのCMOSカメラです。モノクロはきちんとメーカー値とも合うのですが、いまだにカラーの場合でうまく計算できたことがないので、もうここ2年ほど悩み続けています。

Conversion_Factor_ASI290MM_std

同様に横軸が明るさで縦軸がノイズの2乗のN^2になります。測定点が一直線で近似できるのがわかると思います。そのグラフの傾き0.306の逆数1/0.306=3.26がコンバージョンファクターになります。

一眼レフカメラの例は、このページに出てますね。「Details of measurements at each ISO setting」のところからの一連のグラフになります。このようにISO(出力ゲイン)を変えて、順次ISOごとのコンバージョンファクターを測定していきます。コンバージョンファクターが1になるところが「ユニティーゲイン」「ユニティーISO」「ユニティーゲインの時のISO」(言葉だけを知っていても意味がないです、逆に意味をきちんと理解していれば、言葉が多少違っても通じますね)ということになります。

ところが上にも書きましたが、SharpCapではISO(ゲイン)を変えて各コンバージョンファクターを測定することをサボっています。どうせ、コンバージョンファクターはゲインに比例するので、一番低いISOのコンバージョンファクターだけを測って、あとは出力ゲインで割ってやることで、測定回数を劇的に減らしています。これは測定の自動化をするために考えた苦肉策(良く言えば簡単に測定できる方法)といえるでしょう。


読み出しノイズ

これまでのどのグラフでもそうですが、傾きの逆数がコンバージョンファクターになり、信号Sが0の時の切片がその時のISOの読み出しノイズになります。ただし、読み出しノイズに関してはこのページでも書いてあるように
but, for the readout noise it is preferable to measure directly the deviation on a bias image - the result is more precise
と書いてあるように、バイアスファイルから直接測定せよと書いてます。奇しくも、ノイズ会議の時に議論したバイアスと読み出しノイズは同じかというというに対する回答になっていて、やはりバイアス(オフセットとかいう概念と同じ)ファイルで測定されるノイズは読み出しノイズと考えて良さそうです。

読み出しノイズの測定は、カメラにキャップをして最小光量で、時間と共に大きくなるダークノイズを無視するために最小時間で撮影した画像を使います。やはりバイアスの撮影と同じですね。こうやって撮影された画像は読み出しノイズに支配されています。読み出しノイズの直接的な測定についてはこのページを参照してください。
 


測定からいろいろなことがわかる

一連の測定の結果から、非常に重要な幾つかの結論が出ます。例えば、このページで言っている主要なことは
  • Unity gainはISO575のところにある。
  • これは個々のISOについてコンバージョンファクターを測定した結果から来ている。
  • コンバージョンファクターの測定方法は、各ISOで2枚撮影して、その差分から明るさとノイズの関係を評価した。
  • 読み出しノイズはISO6400までは減ってきていて、それ以上のISOでは一定値。なので、ISO6400以上は使っても意味がない。暗い天体はISO6400を使うと(ダイナミックレンジは多少不利になるが、淡いところを出すには)有利
  • 飽和容量は13235ADUと考える(と書いてあが、根拠は不明。14bitだから16348ADUと思ったら、それより小さいので何かかから測定したのか?)。
  • ダイナミックレンジはISO400までは一定で、それ以降は減り始める。なので明るい天体はISO400以下で撮影するのがいい
  • 中間ISOは使うな!例えばISO1000はISO800と同じ読み出しノイズでかつダイナミックレンジは小さい。
ということです。


コンバージョンファクターやユニティーゲインは何の役に立つのか?

答えを一言で言うなら、コンバージョンファクターという(ゲイン依存の)変換係数があるおかげで、ADCの値を読むだけでありとあらゆるものを電子の数で考えることができるようになるので、「単位が揃って便利」だということです。

例えば飽和電子容量(full well)です。本来の飽和電子容量の測定の仕方は、十分にサチレーションを起こすくらい明るい光をカメラに入射し、その時のADCの値の平均値を読み取り、それをコンバージョンファクターで電子の数に変換してやります。コンバージョンファクターが分からなけれが飽和「電子」容量とは言えずに、飽和「ADC」容量とかになってしまいますね。

読み出しノイズの測定もそうです。バイアスファイルは読み出しノイズに支配されています。この時の各素子の明るさのばらつき具合が読み出しノイズになるのですが、当然のことながらこれらの測定は全てADCの出力を見ているので単位は [ADU] で出てきます。こんな時に先に測定したコンバージョンファクターがあると、あーら不思議!なんと電子の数 に変換することができ、普通の読み出しノイズの単位[e-]になります。

逆に言えば、どれだけ画像ファイルからADCでカントされた数を数えることができても、コンバージョンファクターがないと、電子、光子のところまで持っていくことができません。と言うわけでコンバージョンファクターがいかに大切かおわかりいただけましたでしょうか?

では、ユニティーゲインがどうして必要かと言うと、実はそこまで重要な値ではないんですよね。ADU1カウントと測定された電子1個が等しいと言うくらいです。まあ、目安ですね。


電子数と光子数の関係

さらに光子の数sと、センサーで数える電子の数nが、定数で変換できます。この定数をシステム効率ηなどと呼び
\[\eta=\frac{n}{S}\]
と表すことができます。通常はシステム効率は1以下の値をとりますが、ここでは簡単のため1としています。

ポイントはこのシステム効率が内部回路のゲインや積分時間などによらないということです。なので、出てきた電子の数を数えるということは、幾つ光子が入ってきたかが直接わかるため、重宝されるというわけです。

その一方、ADCのカウント数とセンサーで出てくる電子数の関係は内部回路のゲインに依存してしまうため、便利でないのです。

ISOと実ゲインの関係

前回測定していまだに腑に落ちないISOと実ゲインの測定ですが、同様の測定例がどこを探しても見つかりません。ISOとコンバージョンファクターのグラフはすぐに見つかります。これってもしかしたらISOの線形性がほとんどないから出回らないのでしょうか?だとしたら、前回示したグラフは何か失敗しているかと思ったのですが、逆に貴重なのかもしれません。

自分でも各ISOのコンバージョンファクターを測ってみるのが次の課題でしょうか?少なくとも3種類の測定の仕方は考えられるので、それぞれで違いが出るかとかも興味があります。そこで測られたコンバージョンファクターは、やはり実ゲインに比例しているはずなので、もし前回の測定結果のように実ゲインがISOに比例しないなら、どこに矛盾があるか突き止めていくのもまた面白そうです。


おまけ: Unity gainで撮影する意味

unity gainで撮影することには、ほとんど何の意味もないです。これは測定される電子数とADCのカウントの比を表している単なる係数に過ぎません。

たまにunity gainで撮影することが有利だとかいう話を聞くことがありますが、根拠が全くありません。その際によく話されるのが、ADCの1カウント以下で電子(光子でもいい)を測定できないからとかだとかが理由とされますが、そもそも光子を1個だけ数えるのは(量子力学で考えるとあたりまえですが)原理的に不可能です。多数の(単位時間あたりに数にばらつきのある)光子が測定され、統計的に平均値をとると何個の光子が来ていたと言えるだけです。

なので、unity gainに拘らずにISOを決めていいのですが、原理的に考えると最適なISOはDynamic Rangeを損なわない最大のISOということになります。もちろんこれは対象の明るさによってきちんと考えるべきで、明るいもの(昼間の景色や恒星など)がサチるのを避けたいならばISOを下げたほうがいいですし、暗いものを撮影するときはダイナミックレンジを犠牲にしてでもISOをあげた方がいい時があります。


まとめ

コンバージョンファクターについて、できうる限り簡単に書いてみました。みなさん理解できましたでしょうか?わかりにくいところとかありましたら、コメント欄にでもお書きください。

もう少し6D測定を続けてみたいと思います。他の結果と矛盾がないのか、それとも何か間違っているのか?どのような設定で撮影すればいいかの根拠になっていくので、これはこれでかなり楽しいです。



とうとう念願のEOS 6Dのユニティーゲイン(unity gain、電子とADCの1カウントが等しくなるゲイン)を測定してみました。といってもまだ思うところもあるので、暫定的な結果です。


これまでの経緯

使ったのはSharpCapで、昨年9月ころの3.3βから一眼レフカメラをサポートし出したため、もしかしたらLive view機能でシャッタを切り続ければ、センサー解析機能を使って一眼レフカメラのセンサーも解析できるのではと思ったからです。



ShapCapを使ったセンサーの解析手法についてはASI294MCなどを測定していて、メーカー値とかなり一致することがわかっています。





6Dでの測定

さて、実際に6DをSharpCapに繋いで、「センサー解析」を使用してみましょう。使ったSharpCapは2021/3/10リリースの最新の4.0.7493.0(BETA)の64bit版です。

ところで、なんでわざわざカギ括弧付きでセンサー解析と書いたかというと、メニューとカメラ制御の部分の日本語化に貢献しているからです。いのさんと智さんも貢献してくれました。特にいのさんは私の拙い訳をかなりまともな用語に直してくれました。

さて、まずセンサー解析を立ち上げますが、やはりどうも「ライブビュー」モードにしないとそもそも機能しないみたいです。逆にいえばライブビューモードにさえしておけば、あとはほとんどCMOSカメラと同じ操作になりました。

まず大事なのは光源の明るさ設定。目標はBrightnesのとこの50%付近に鋭いピークが立つこと。私は以前と同様にiPadのColor Screenというアプリを使い、Hue0、Saturation0、Brihgtness 32となるようにして、6Dのレンズを外しそのままiPadの上にセンサー面が向くように置きました。エリア選択がありますが、選択範囲内で周辺減光など光のスロープがあるとノイズが必要以上に大きく出てしまうので、センター付近の比較的狭いエリアを選びます。円状のレチクルを出して、その一番小さい円に内接するように正方形のエリアを決めました。あとは初期の光の量がまずいと怒られるので、ISOを100、露出時間を250msとしたら解析スタートの許可が出たので、そのまま進めます。

IMG_1980
セットアップの様子。

あとはひたすら待つだけです。CMOSカメラと違い、1フレームづつ撮っていくので時間とコマ数がかかります。終了まで約1時間ちょっと、1000回弱のシャッターを切りました。SharpCapからASCOMを通じて6Dの露出時間とISOを随時切り替えてシャターを切ります。ただしSDカードに記録はしないため、バッテリーは1000枚撮影したあともまだフルゲージ残ってました。

IMG_1972
下にフレーム数と時間が出ています。

今回測定したゲインはISO100からISO500までの8段階でした。それぞれのゲインで、センサーの読み取りから、暗すぎたりサチったりしないように適当に露出時間にフィードバックして適した露出時間を決めるため、それだけで何度もシャッタを切るので、どうしてもシャッター回数が多くなってしまいます。

一眼レフカメラのシャッター回数は寿命に繋がるので、無駄な機械シャッターを切らないように少なくとも何度かCMOSカメラで練習することをお勧めします。

以前のバージョンの測定の時には、この適した露出時間がなかなか決まらなくて長くしたり短くしたりを永遠と繰り返すバグなどもありましたが、今回はそのようなことはなかったです。ただし、測定中に露出時間も同じでISOも同じなのに撮影した画面に出てくる明るさがあからさまに変わって、安定しないような時がありました。原因はわかりませんが、ここは少し結果に対して不安要素となっています。

途中ダークノイズの測定のためにキャップをしたり、終わったら外したりしますが、それらは指示に従えばいいでしょう。 


測定結果

結果を見てみます。

IMG_1977

Gain Valuee/ADURead Noise (e)Full Well (e)Relative GainRel. Gain (db)Dynamic Range (Stops)
1005.6927.18931441.000.0011.74
1254.4526.69730541.282.1111.42
1603.2912.66730541.734.7412.06
2002.4911.87407922.287.1711.75
2501.9011.393107239.5411.41
3201.365.77223064.1812.4111.92
4000.954.99155525.9915.5511.60
5000.744.87121217.6817.1111.28

グラフ化しておきます。一番下のグラフは読み出しノイズをADUで表した場合です。ISO300くらいまではゲインが上がると共にe-単位での読み出しノイズが小さくなっていくのでほぼ一定で、ISO300を超えるとADUで見て読み出しノイズが上がってきます。これはISO300を超えると実際の画像で見て読み出しノイズが大きく見えてくるということを示しています。

6D_gain_graph_cut


これだけみていると、unity gain (unity ISO)は e/ADUが1になるところなので400を切るところ程度と読めます。ところがこの値は少し疑問が残ります。

Read noiseについては3段階に分かれているようなので、ここから3種のアナログゲインがあるのではとかの推測ができます。さらに、細かいゲインについてはデジタルゲインの可能性が高いと言えます。


考察

まずはこちらのページを見てください。


6Dのセンサーについて測定しているページです。このページではunity gainはISO575であると言っています。今回自分で測定したISO400弱とというのとは1.5倍くらいのズレがあります。

わかりやすくするために上記ページの表のISO800までをお借りします。

ISOgain[e-/ADU]read noise[e-]DR[dB]
505.641327.4868.7
1005.691127.5868.7
2002.773213.6668.6
4001.41787.5367.9
8000.72914.4566.7

Read noiseについてはISO100、200、400、800にアナログゲインが入っていると考えられるので、今回測定した結果と矛盾ないと考えられます。

gain[e-/ADU]については少なくとも今回測定した値の中でISO100のところはgainもread noiseも上の測定結果とよく合っていると言っていいと思います。ところがそれ以外のISOのところは全て1.5倍くらいずれています。これはどういうことなのでしょうか?

この違いは、今回SharpCapが簡易的な測定をしていることに起因します。先に示して表の中で、実際に測定しているところを赤くしてみます。

Gain Valuee/ADURead Noise (e)Full Well (e)Relative GainRel. Gain (db)Dynamic Range (Stops)
1005.6927.18931441.000.0011.74
1254.4626.69730541.282.1111.42
1603.2912.66539621.734.7412.06
2002.4911.87407922.287.1711.75
2501.9011.393107239.5411.41
3201.365.77223064.1812.4111.92
4000.954.99155525.9915.5511.60
5000.744.87121217.6817.1111.28

この赤いところ以外は実測ではなく、実測した値から計算しているに過ぎません。例えば
  • e/ADUのISO100の5.69以外のところは、5.69を単にRelative Gainで割った値に過ぎません。
  • Full Wellも同様で、ISO100のところの93144を単にRelative Gainで割った値に過ぎません。
  • Dynmic Rangeは計算値のFull Wellを実測のRead Noiseで割ってbitで表しただけです。
こうやってみると、Dynmic Rangeが今回測定した範囲の中であまり変わらないのも理解できます。なぜなら、Read Noiseがアナログゲイン(ISO)に比例してよく下がってくれている範囲内だからです。今回の測定範囲外は上記ページを見てもらえばよくわかります。ある一定値以上のISOでは、これ以上いくらアナログゲイン(ISO)をあげても、Read Noiseが他の要因である一定値に制限されてしまう一方、Full Wellは下がっていくために、ダイナミックレンジが小さくなっていきます。

また、Dynamic Rangeで考えたら、最大値とほとんど変わらないISO1000位までまではISOを上げたほうが得。Dynamic Rangeの落ちを1bitまで(半分になるということ)許容するとしたら、ISO1600までは許容範囲で、ISO3200だとそれよりほんの少し損をするといったところでしょうか。なので私がもし使うならISO800か1600、もう少し明るさが欲しい場合はぎりぎりISO3200ということにするのかと思います。


今回の測定の問題点

さてここで、今回実測したRelative Gainのところに注目します。通常はISO100を基準にISO200なら2倍、ISO400なら4倍になるはずですが、結果はISO400で2.3倍、ISO800で6倍でどうも1/1.5倍程度小さく測定されてしまっているようです。しかも線形性もあまりないという冴えない結果です。先に、測定中に撮影した明るさが一定にならないと書きましたが、これが悪さをしている可能性があります。もしISOと実際のゲインが理論値で一致しているなら(ISO200なら2倍、ISO400なら4倍とかいうこと)unity gainはISO569になるはずで、上記ページの結果ともほぼ一致します。

SharpCapでの測定方法は、CMOSカメラが前提のためにシャッター回数を気にしなくてため、何枚も画像を撮り、それらの同じエリアを使うことで、時系列でずれたようなデータを使い解析しています。一方、1枚の暗いところから明るいところまで含まれるような画像を撮影し、その画像の中で空間的にずれたようなデータを使うことでも同様の解析ができます。

というか、普通は後者の方が素直なやり方で、SharpCapは測定を自動化するために時系列のデータを使っているというわけです。最初SharpCapのやり方を見た時に「上手いやり方だなあ」と思いましたが、測定してみると簡易的な方法であることはすぐに認識できて、しかも今回一眼レフカメラのシャッター回数のことを考えると、やはり1枚の画像で解析した方がいい気がしてきました。

明るさのところをもう少し改良して、もう一度くらいなら測定したいと思います。シャッターの寿命が10万回としたら、1回の測定でシャッター寿命の約1%を使ってしまう計算なので、SharpCapで測定するのは最小限に抑えておいたほうがいいでしょう。

むしろISO100の結果だけは正しく測定していて、結果も矛盾なく正しく得られていると思われるので、あとは自分で別途ゲインを測定したほうがマシそうです。


ISOと実ゲインを実測

というわけで、ISOと実際に撮影できる明るさを実測してみたいと思います。

測定方法はSharpCapで測った時と同様に、iPadのアプリColor Screenを使い、そこにカメラを載せて、ISOを変えて撮影します。ただし、SharpCapでの測定がばらついたのでその反省を生かし、外光の影響ができるだけないのようにしました。まず、Color Screenの明るさを32から128の4倍にします。さらに、Color Screenの上に薄い紙を一枚敷きiPadの表面の反射の影響をなくすようにします。先の測定ではレンズなしのセンサー面を暴露しての測定でしたが、これもレンズをつけてレンズの先の光だけがセンサーに入るようにしました。

露光時間を1/50秒に固定して、ISOを800から100まで下げて撮影していきます。というのはISO800でサチらないように気をつけるためです。

測定は、撮影した各画像の中心の100x100ピクセルの平均の明るさ。RGB個別に取り出してます。中心を選ぶ理由は、縁のほうに行くと周辺減光が影響してくるからです。また、ADCの値に何らかのオフセットが加わっている可能性があるために、レンズに蓋をして明るさを測りましたが、レンズを開けた時に対して0.5%程と測定にほぼ影響はないので今回は無視しました。

これだけやったのですがやはり結果は変わらず、ISOと実際のゲインは全然合いませんでした。結果を示します。ISO100の時のゲインを1としています。
6D_ISO_gain
普通はISO800ならISO100の8倍の明るさになるはずです。グラフの点線に近くなるはずです。ところが実測は4-5倍程度しかないどころか、線形性(測定点を結んだ線が真っ直ぐになること)さえもありません。

測定はかなりしっかりやったつもりです。いったい何がおかしいのでしょうか?一つの可能性は、ヒストグラムのピークが一定の場所になるように露光時間を変えるなど調整しながら、その露光時間ぶんを補正してISOを変えて測定するとかでしょうか?ちょっといろいろ不明で、そろそろ力尽きたのでここは次の課題とします。


今回の結論

はっきりとした結論は出ませんでしたが、まとめます。
  • SharpCapのセンサー解析機能を使うことで、EOS 6DのISO100についてのコンバージョンファクター(Gain)はきちんと測定できたようです。
  • ですが、それ以外のところは1.5倍ほどずれていると考えられます。
  • その原因は、ISOを変えることに対して明るさが期待通りにならないことから来ていると思われます。
  • なので、unity gainに関しては保留とします。
  • 他の場所ではunity gainはISO600を切るくらいと、複数確認できるのと、ISOとゲインの関係さえしっかりしたら今回の測定でもそのくらいになるので、おそらくunity gainはISO600を切るというのが正しいと思われます。
今ふと思ったのですが、もしかしたら持っている6Dのゲイン設定がおかしくなっていて、本当にISOとずれているのかもしれません。こうなったら修理コースなので、もう少しいろいろ試してから結論を出したいと思います。


あぷらなーとさんはじめ何人かの方がMatlabを購入し画像処理に活用し始めましたようです。



私もMatlabは学生の頃から使っていて、今調べてみたら一番古いファイルは1998年でした。なので20年以上使っていることになりますが、 これまで画像処理に使ったことはありませんでした。

あぷらなーとさんが指摘しているように、確かにMatlabは行列を容易に扱えるので、2次元が基本の画像処理は向いているはずです。最近過去のコードをpython化したりしてるのですが、pythonは2次元だと必ず繰り返し処理が入るので、コードが冗長になりがちです。私も少し試してみたのですが、どうやら画像処理に関してはMatlabの方がコードがかなりシンプルになり見通しが良くなるので、有利というのは正しそうです。とうわけで、遅ればせながら私もMatlab画像処理に参画します。

Matlabを使ってまず手始めに試したのが、昨年3月依頼謎が解けなくて止まってしまっているASI294MC Proのノイズ解析の続きです。

 


最初にごめんなさいと謝っておきます。今回の記事はかなり細かくてめんどくさい話で、ほとんどの人には役に立たないです。しかもうまくいかなかったという結果なので、ほんとに興味のある方だけ読んでみてください。

Matlabで画像処理を始めてみたいという人がいたら、もしかしたらコードが参考になるかもしれません。


以前のおさらい 

ざっとおさらいしてみます。ZWO社のASIシリーズは性能を示すデータを公開してくれています。例えば広いセンサー面積を利用した電視観望や、撮影にも十分に使うことができるASI294MC Proのページを見ていくと、下の方にいくつかグラフが出てきます。このグラフ、一見分かりそうで分かりにくいところもあるので、以前その読み方を解説しました。



上のページでも解説している通り、SharpCapを使うとZWOのページにあるようなデータを測定することができます。冷却バージョンのProも



で実際に測定した記事を書いています。特にGain[e/ADU]の値はコンバージョンファクターとも呼ばれて、設定gainが0の時のコンバージョンファクターはユニティーゲインなどにも直結するような非常に重要な値になります。少しコツは必要ですが、SharpCapのスクリプトでCMOSカメラの特性を実測すると、ZWOにあるような値に非常に近いものを測定することができます。2つのツールで同様のデータが取れているということはかなり信頼が置けるはずです。さらにもう一歩進めて、このようなツールでのスクリプトでなく、実際に自分でノイズを撮影して解析してみようと試したのが、昨年3月までの一連の記事になります。

ところが、実際に撮影して解析してみるとどうしてもZWOのデータやSharpCapでの解析結果と合いません。

SharpCapで撮影した時の撮影条件は限りなく再現させたと思っています。具体的には撮影対象の明るさ、カメラのゲイン、露光時間です。明るさは同等のものが撮れているので、撮影に関しては問題ないと思います。問題はノイズです。明るさとノイズの関係からコンバージョンファクターを求めるのですが、撮影された画像の明るさのばらつきが、想定していたものよりかなり大きいと出てしまいます。

具体的には標準偏差を用いるとノイズが大きすぎると出てしまい、苦肉の策で(ノイズが小さいとでる)平均偏差を使うと、大体ZWOとSharpCapの測定と一致するという結果でした。

ヒントになるかと思いモノクロのASI290MMで測定したら、標準偏差で計算してもZWOやSharpCapと一致した結果を得ることができました。ということはやはりカラーの場合も平均偏差などを使わずに標準偏差を用いるのが正しいのではと推測することができます。

そんな折、あぷらなーとさんがRGGBを一素子づつ解析したCFA(Color Filtr Array)で評価した場合と、RGBにdebayerした場合では、debayerの方が標準偏差で考えて0.75倍とか0.79倍程度に小さくなることを示してくれました。



それでもdebayerで計算してもまだZWOやSharpCapでのノイズの少なさには辿りつきません。

いろいろやっていて結局行き着いたところはこの画面で、

IMG_6436

小さいと思われるRGBの標準偏差よりも、Lの標準偏差の方がなぜかさらに小さいことがSharpCapの一枚の撮影で出てしまっていることです。いまだにSharpCapが内部でどうやって計算しているのかよくわかりません。このLの標準偏差を仮定するとノイズはかなりZWOもしくはSharpCapのスクリプトで測定した結果に一致します。言い換えると、これくらい小さい標準偏差くらいのばらつきになっていないと結果は一致しないということです。


Matlabでの画像処理の実際

やっと前振りが終わりましたが、今回Matlabで以前のpythonコードを書き直すかたらわら、どうしたら一致した結果を出せるか、なぜSharpCapのLがノイズが小さく出ているかを念頭に置きながらいろいろ試してみました。

Matlabを初めて使う人が一番面食らうのは、配列の表記法かと思います。これが独特なのですが、この独特な表記によりコードがシンプルになるので避けるわけにもいきません。私が書くよりももっといい説明がたくさんあるかと思いますが、今回の画像処理に必要な最低限だけ書いておきます。

まず2次元の配列、Matlabだと行列といっていますが、例えば2行x3列の配列Aを

>> A=[1 2 3;4 5 6;]

A =
     1     2     3
     4     5     6

と作ります。例えば第2行,第1列成分を見たい場合には

>>A(2,1)

と打つだけです。答えは

ans = 4

と出ます。これは至って普通です。独特なのは:(コロン)の使い方で、例えばAの1行目すべてを見たい場合は

>>A(1,:)

と打ちます。答えは 1 2 3となります。:はどこからどこまでを意味し、

>>A(1,1:2)

だと

ans =     1     2

となります。:(コロン)だけだとその次元の最初から最後まで全部という意味です。これがMatlabで絶対理解しておかなければならない表記法の一つです。

今回は画像を扱うのに4次元の配列を使っています。1、2次にy座標とx座標。左上からyは向かって下に、xは右に移動していきます。3次目はRGBやCFAのインデックス。RGBなら3つ、CFAなら4つとっています。 4次目は何枚目の画像かを表しています。今回8枚の画像を使っていますが、その何枚目かを表します。あと成分を表す数字は1から始まることにも注意です。0からではありません。 

例えば、AにRGBで別れた画像が8枚入っているとすると、5枚目のG画像を表す時は

A(:, :, 2, 5)

と表します。:が2つあるのはy、x座標それぞれ全部を表しています。"2"はRGBの2番目と言う意味です。Bなら3ですね。"5"は5枚目と言うことです。

例えばサンプルコードの最初の方にある

 Raw(:, :, i) = fitsread(ファイル名);

はi番目の画像にfitsファイルを読み込んで、1次目にy座標のデータ、2次目にx座標のデータを全部入れていきなさいと言うのをわずか1行で表記したものです。

これらの表式は慣れるしかないのですが、ここらへんを感覚的に理解できるとコードもすらすら読めて、シンプルなコードを書くことができるようになるのかと思います。


モノクロセンサーASI290MMの場合

とりあえずMatlabでやってみた画像処理を説明していきます。まずはモノクロのASI290MMの結果です。

ASI290MM
このグラフを出したMatlabのソースコードをアップロードしておきます。使用した画像も(サイズを切り詰めるために中心部のみ切り取ったものですが)一緒に入れておきましたので、すぐに走らせることができると思います。もしMatlabを持っていて興味があったら走らせてみてください。

ASI290MM.zip 


結果は前述したとおり、平均偏差ではなく一般的な標準偏差を使っています。メーカー値はコンバージョンファクターが3.6程度、ユニティーゲインが110程度ですので、グラフ内の数値と比較してみるとそこそこ一致していることがわかります。なので、カラーでも標準偏差を用いる方向を探りたいと思います。

また、以前pythonで書いたコードと比較して同等の結果を出すことできています。

Conversion_Factor_ASI290MM_std

Maltabでもこれまでと同様のことができることがわかり、かつ遥かにシンプルに書けることがわかりました。


カラー版ASI294MC Proの場合: CFAで解析 

さて、Matlabでの画像処理にも慣れたところで問題のカラーのASI294MC Proの解析です。

結論から言うと、結局今回も謎は解けずじまいでした。

まずはシンプルなCFAの方から。RAW画像からCFAへ分離するところはビルトインの関数とかはないようなので自分で書いてみました。以前は平均偏差(mad)で無理やり答えを合わせましたが、今回ASI290MMの結果なども含めていろいろ考えた末に、結局諦めて標準偏差(std)を使うことにしました。

ASI294MCPro_CFA

各測定点はそのまま解析した結果を表示していますが、フィッティングは無理やり合うように補正してます。考え方としては、謎な部分をunknown factorとしてフィッティングするところで無理やり補正するというやりかたです。今回のCFAの場合、unknown factorはノイズ(標準偏差の2乗)を2分の1とすることでZWO、SharpCapの結果と一致することがわかりました。

ちなみにメーカー値はコンバージョンファクターが3.9程度、ユニティーゲインが117程度です。繰り返しますが、CFAで測定されたノイズを無理やり半分にすることでメーカー値に近くなります。

先に説明したように、SharpCapでLのノイズが小さく出ている理由がわからなかったので、もうこのように無理やり合わせてしまうしかなかったという苦肉の策です。これはあぷらなーとさんが示してくれた、CFAとRGBで標準偏差で0.75倍違うというのを考慮してもまだ少し足りません。

まあ、このやり方に言いたいことはたくさんあるかと思いますが、とりあえず先に進むことにします。


カラー版ASI294MC Proの場合: RGBで解析

次はRGBでの解析です。

最初はMatabビルトインのdemozaicという関数を使ったのですが、これだと中で何か変に処理しているのかノイズが話にならないくらい大きくなりすぎました。仕方ないので自分でRGBに落とすコードを書いています。そこそこまともそうな値は出たのですが、ただしこのコードまだ未完成です。なぜかというと、センサーアレイはRGGBなのでG(グリーン)が二つあるのですが、その応答が少し違うことに起因します。Gを2つ使って求めると強度の山が2つでできてしまい、ばらつきが大きくなりすぎるからです。そのため今回は2つあるG素子のうち1つだけを使うことにしました。

ASI294MCPro_RGB

RGBになりあぷらなーとさんが示してくれた通り、CFAの時よりもばらつきは少なくなることがわかりました。それでもまだノイズが大きすぎるのでCFAの時と同様にunknown factorとしてフィッティングするところに無理やり入れました。RGBの場合、unknown factorはノイズ(標準偏差の2乗)をルート2で割ってやることでZWO、SharpCapの結果とかなり一致することがわかりました。


考察

ASI294MC Proで使ったファイルもアップロードしておきます。よかったら参考にしてください。

ASI294MCPro.zip 

今回はまだ苦肉の策で、無理やり答えを合わせるように測定結果を割っているだけなので、考察と呼べるようなものにはなりません。CFAの場合出てきたノイズを半分に、RGBの場合ルート2でわってやるとSharpCapやZWOの結果とかなり近い結果になりましたが、その理由は全然わかっていません。

前回のpythonコードを使った時はCFAやRGBの変換はPixInsightを使いました。でも今回はRAWファイルから直接計算してCFAもRGBも取り出しています。ここら辺はMatlabに移ったことでやりやすくなったところだと思います。このように今回はかなり噛み砕いてブラックボックスがないように進めてきたつもりです。撮影した画像はSharpCapのスクリプトで撮影したものと同等と仮定していいかと思います。これでも結果が一致しないということは、
  • ZWOもSharpCapも何か特殊なことをやっているか
  • もしくはまだ私が馬鹿なミスを犯しているか
です。とにかく怪しいのがZWOのLのノイズが小さく出過ぎていること。この謎が解けない限り進展がない気がします。


まとめ

週末を結構な時間費やしましたが、コンバージョンファクターに関しては結局昨年から進展がなかったので疲れてしまいました。それでも得られた成果としては、
  • Matlabでの画像解析の手法に慣れることができた。
  • あぷらなーとさんが出してくれた、debayerの方がCFAよりもノイズが0.8倍程度と小さくなることを確かめることができた。
ことくらいでしょうか。でも全然だめですね。これ以上は時間の無駄になりそうなので、一旦ここで打ち切ることにします。

それでもMatlabが使いやすいことは分かったので、今後の画像解析はMatlabで進められそうな目処はついたと思います。

元々この解析はダークノイズ解析につなげたいとか、EOS 6Dのユニティーゲインを測定したいとかで始めたのですが、なかなか進みません。いつ実現できるのかまだまだ未定ですが、諦めることはしたくないのでまたいつか再開する予定です。

このページのトップヘ