ほしぞloveログ

天体観測始めました。

タグ:RGB

あぷらなーとさんはじめ何人かの方が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のユニティーゲインを測定したいとかで始めたのですが、なかなか進みません。いつ実現できるのかまだまだ未定ですが、諦めることはしたくないのでまたいつか再開する予定です。

連休ですが、富山は相変わらず天気はよくありません。こんな時はたまっていた課題を片付けます。


CMOSセンサーを顕微鏡で見るとどうなる?


事の発端は、小海の星と自然のフェスタでスタベのアルバイトのK君がASI294のセンサー面を実体顕微鏡で見てみたいと言った事です。ASI294は電視観望用に持っていますが、うちにはたまたま実体顕微鏡もあります。下の子Sukeの自由研究のために2017年の原村星まつりで買ったものです。これならなんとかなりそうです。

ところがシベットさんからコメントで実体顕微鏡では倍率が低すぎるのではという指摘がありました。なぜかうちには倍率の高い生物顕微鏡もあります。これも同じく下のSukeの自由研究のために2016年の原村星まつりで購入しています。というわけで、必要なものはそろっているので実際に見ることにしました。果たして何が見えるのか?

IMG_9068


まずは実体顕微鏡

手持ちの実体顕微鏡は対物レンズが4倍、接眼レンズが10倍で、40倍のものです。ここにASI294を置いてやればいいのですが、カメラ本体が大きすぎでレンズに近づき過ぎてしまいピントが出ません。仕方ないので、今回はもっと奥行きの短いASI224MCで試すことにします。K君、ごめんなさい。

ASI224MCのセンサーサイズは4.8×3.6mm、解像度は1304×976で一素子のサイズは3.75×3.75ミクロンとのことです。

IMG_9091

センサーは思ったより小さく、周りに金色の線がたくさんつながっています。これを40倍の実体顕微鏡で見てみました。とりあえずiPhoneのカメラで撮ってみます。

IMG_9036

よく見えていますが、センサー面はのっぺりしているだけでやはり全然倍率が足りないことがわかります。

でも実体顕微鏡は簡単には倍率を変えることができません。いろいろ考えて、手持ちのアイピースを利用することにしてみました。接眼レンズを片方外します。31.7mmのアメリカンサイズは径が大き過ぎたので、昔の1インチタイプのものを使いました。手持ちは20mmと12.5mmと9mmの3つ。

IMG_9071

最初に20mmのもので見てみます。1インチサイズだと今度は径が小さすぎるので、とりあえず固定せずに穴に入れるだけです。

IMG_9039

もともとついていた10倍のレンズより少しだけ大きく見えます。

IMG_9037

10倍レンズは焦点距離25mmとか30mmくらいなのかと思います。この時点で50倍とか60倍相当かと思います。当然これだけだとまだ倍率は不十分なので、さらに倍率を上げるために12.5mmレンズを取り付けます。

IMG_9041

これも穴に入れただけで固定していませんが、まあとりあえずそこそこ安定して見えるものです。像はというと、

IMG_9029

接続線は大きく見えてきましたが、センサー面はほとんど変化なしです。最後9mmを試します。

IMG_9043

IMG_9061

ここでiPadのLEDで光を当てたこともあり、基板のパターンが見えてきました。しかもピントを調整してセンサーの端をよく見ると目盛りのようなものが見えます。

IMG_9064

小さな目盛り10個につき大きな目盛りがあり、大きな目盛りを数えると長辺で14個、短辺で10個程度あります。センサーサイズが4.8×3.6mmなので、大きな目盛り2個で1mmはないくらいということがわかります。上2枚の写真のセンサーの枠の太さは同じなので、遡っていくと全体サイズまで直接比較することができると思います。

これで倍率は手持ちの1インチアイピースでは限界です。最終倍率は120倍程度かと思われますが、それでもセンサー面の構造は全く見えてきません。やはりもっと倍率の高い生物顕微鏡が必要なようです。


生物顕微鏡なら見えるか?

IMG_9094


さて、生物顕微鏡に移ったのですが、ここで問題が発生しました。ASI224MCでも奥行きがあり過ぎて、ステージを一番下に下げても対物レンズの間に入らないのです。

IMG_9093


最初はセンサーがついている基板を外そうかと思いました。

IMG_9092

ですが、基板の右側についているコネクタをうまく外すことがどうしてもできなくて、泣く泣く諦めることにしました。ちなみに左上についている赤い検品シールは触るとすぐに崩れるタイプで、ネジを緩めるとシールが簡単にとれてしまいます。改造したかどうかのチェックを兼ねている様です。なのでこのカメラはもうメーカーのサポートは受けることができないと思っておいたほうがよさそうです。もし個人で分解する際は気をつけてください。

結局今回は仕方ないので、下の写真の様に顕微鏡のステージを取り外すことにしました。

IMG_9076

ステージに置くことができないので、CMOSカメラを手で持ちながら見ることになります。バランスが悪いのですが、何度か撮影すればいい位置とピントで写ることもあり、何とかなりそうです。

倍率は対物レンズが4倍、10倍、100倍。接眼レンズが10倍と15倍のものがあります。

とりあえず対物4倍と接眼15倍で目で見てみます。明かりがもっと欲しかったので、先日中身を取り替えたパワータンクの強力LEDで照らすことにしました。この状態で目で見てみると、手でカメラを押さえながらですが、何とか見ることができます。

ここで秘密兵器登場。今回の撮影のために、SVBONYの顕微鏡の差込口に挿すことができるアダプターを買っておきました。

IMG_9089

これはT2ネジが切ってあるので直接ZWOのASIカメラに取り付けることができます。

IMG_9084

これを生物顕微鏡の接眼レンズの代わりに取り付け、直焦点撮影をします。

最初はASI294MC Proを撮影用のカメラとして使いました。まずは最低倍率の4倍の対物レンズで撮影します。撮影はMacのASICAPを使いました。

2020-01-12-0626_0-CapObj_0001

すでに先ほどより大きく写り、遥かに精細に写っています。だんだん楽しくなってきました。

対物レンズの倍率を4倍のものから10倍のものに上げます。

2020-01-12-0607_4-CapObj_0000

おお!とうとうセンサー面の構造が見えてきました。

もっと倍率をあげたいのですが、100倍の対物レンズはセンサー面に相当近づけなければならなく、保護カバーを外さなければピントが合わなさそうです。さらに、かなり暗くなることがわかっているのであまりやりたくありません。いろいろ考えて、ASI294MCよりも一素子のサイズが小さくて分解能の高いASI178MCを使ってみることにしました。その結果がこちらです。

もうセンサー面の構造もはっきりと見えます。

2020-01-12-0634_3-CapObj_0000

でも構造は見えますが、RGBフィルターとかの影響がわかりません。目盛りがついているところから少し中に進むと紫のエリアがなくなる境目があるのですが、わかるのはこれくらいです。この頃になるとカメラの固定方法の工夫でかなりピントも合わせやすくなってきました。下はピントや撮影カメラのゲインを相当合わせ込んだ場合です。

2020-01-12-0636_7-CapObj_0000

かなりはっきり見えて、センサー10素子で目盛りが一つ進むこともわかります。なので目盛りはセンサーの数を表していたんですね。でも、これでも全部同じ素子のように見えてしまっています。

どうやっても進展がないので、ここでかなり悩みました。最後にとった方法が、わざとカメラを回転させて素子を45度傾けてみること。理由ははっきりとわかりませんが、これでやっとうまく見ることができました。最終結果です。

2020-01-12-0702_0-CapObj_0000

見事にRGGBフィルターの配列を見ることができます。これを見て改めて思うのですが、これだけの微細構造を作る技術は見事なものです。これが民生レベルで安価に購入できるというのはなんと幸せな世の中なのでしょう。

最後に、その時の撮影風景です。

IMG_9083


考察

とりあえずセンサー面の構造とRGGBフィルターの存在を確認することができました。対物レンズの倍率を上げ、十分な灯りを用意することができれば、もっと微細な構造を見ることができるかもしれません。

RGGBフィルターがかかっている部分と、かかっていない部分の境目もはっきりとわかります。RGGBフィルターは全面にかかっているわけではなく、センサーの端の方はフィルターはないことがわかります。


問題はここから何を引き出すかです。K君と話した時、サッポロポテト現象を解明したというようなことを言っていた気がします。サッポロポテト現象はあぷらなーとさんも困っているようです



サッポロポテト現象は各素子についているマイクロレンズ効果が原因のようですが、はっきりとした解決策はあまりないようです。今回はやっとCMOSセンサーの素子の構造が見えたくらいなので、まだゴールまでは程遠いと思います。

今回の結果をどう活かすのかが今後の課題でしょう。


まとめ

こんなふうに工夫でどんどん見えてくるような実験は超面白くて大好きです。

とりあえず今回分かったことは、
  • 実体顕微鏡ではCMOSセンサーの構造を見るには倍率が低すぎる。
  • 生物顕微鏡を使うことでCMOSセンサー面の構造、RGBフィルターの様子を見ることができる。
  • センサー1素子そのものを十分な解像度で見るには至らなかったが、まだ拡大率を上げる手は残されているので、1素子をもう少しはっきり見ることもできるかもしれない。

課題としては
  • センサーを置く架台をしっかりしたほうがいい。
  • センサー面についている保護ガラスを外して、高倍率で撮影してみる。
  • 動画撮影でスタックするとさらに解像度が上がるかもしれない。
  • ASI294MCのセンサー面を見てみる。
ことくらいでしょうか。あー、今回も面白かった。


K君、遅くなってすみませんでした。やっとなんか見えるくらいにはなってきました。

K君、次は何を見てみたいですか?


このページのトップヘ