あぷらなーとさんはじめ何人かの方が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でのノイズの少なさには辿りつきません。
いろいろやっていて結局行き着いたところはこの画面で、
小さいと思われるRGBの標準偏差よりも、Lの標準偏差の方がなぜかさらに小さいことがSharpCapの一枚の撮影で出てしまっていることです。いまだにSharpCapが内部でどうやって計算しているのかよくわかりません。このLの標準偏差を仮定するとノイズはかなりZWOもしくはSharpCapのスクリプトで測定した結果に一致します。言い換えると、これくらい小さい標準偏差くらいのばらつきになっていないと結果は一致しないということです。
やっと前振りが終わりましたが、今回Matlabで以前のpythonコードを書き直すかたらわら、どうしたら一致した結果を出せるか、なぜSharpCapのLがノイズが小さく出ているかを念頭に置きながらいろいろ試してみました。
Matlabを初めて使う人が一番面食らうのは、配列の表記法かと思います。これが独特なのですが、この独特な表記によりコードがシンプルになるので避けるわけにもいきません。私が書くよりももっといい説明がたくさんあるかと思いますが、今回の画像処理に必要な最低限だけ書いておきます。
まず2次元の配列、Matlabだと行列といっていますが、例えば2行x3列の配列Aを
>> A=[1 2 3;4 5 6;]
と作ります。例えば第2行,第1列成分を見たい場合には
>>A(2,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枚目と言うことです。
例えばサンプルコードの最初の方にある
はi番目の画像にfitsファイルを読み込んで、1次目にy座標のデータ、2次目にx座標のデータを全部入れていきなさいと言うのをわずか1行で表記したものです。
これらの表式は慣れるしかないのですが、ここらへんを感覚的に理解できるとコードもすらすら読めて、シンプルなコードを書くことができるようになるのかと思います。
とりあえずMatlabでやってみた画像処理を説明していきます。まずはモノクロのASI290MMの結果です。
このグラフを出したMatlabのソースコードをアップロードしておきます。使用した画像も（サイズを切り詰めるために中心部のみ切り取ったものですが）一緒に入れておきましたので、すぐに走らせることができると思います。もしMatlabを持っていて興味があったら走らせてみてください。
ASI290MM.zip
メーカー値はコンバージョンファクターが3.6程度、ユニティーゲインが110程度ですので、グラフ内の数値と比較してみるとそこそこ一致していることがわかります。
また、以前pythonで書いたコードと比較して同等の結果を出すことできています。
Maltabでもこれまでと同様のことができることがわかり、かつ遥かにシンプルに書けることがわかりました。
さて、Matlabでの画像処理にも慣れたところで問題のカラーのASI294MC Proの解析です。
結論から言うと、結局今回も謎は解けずじまいでした。
まずはシンプルなCFAの方から。RAW画像からCFAへ分離するところはビルトインの関数とかはないようなので自分で書いてみました。以前は平均偏差(mad)で無理やり答えを合わせましたが、今回ASI290MMの結果なども含めていろいろ考えた末に、結局諦めて標準偏差(std)を使うことにしました。
各測定点はそのまま解析した結果を表示していますが、フィッティングは無理やり合うように補正してます。考え方としては、謎な部分をunknown factorとしてフィッティングするところで無理やり補正するというやりかたです。今回のCFAの場合、unknown factorはノイズ(標準偏差の２乗)を2分の1とすることでZWO、SharpCapの結果と一致することがわかりました。
ちなみにメーカー値はコンバージョンファクターが3.9程度、ユニティーゲインが117程度です。繰り返しますが、CFAで測定されたノイズを無理やり半分にすることでメーカー値に近くなります。
先に説明したように、SharpCapでLのノイズが小さく出ている理由がわからなかったので、もうこのように無理やり合わせてしまうしかなかったという苦肉の策です。これはあぷらなーとさんが示してくれた、CFAとRGBで標準偏差で0.75倍違うというのを考慮してもまだ少し足りません。
まあ、このやり方に言いたいことはたくさんあるかと思いますが、とりあえず先に進むことにします。
次はRGBでの解析です。
最初はMatabビルトインのdemozaicという関数を使ったのですが、これだと中で何か変に処理しているのかノイズが話にならないくらい大きくなりすぎました。仕方ないので自分でRGBに落とすコードを書いています。そこそこまともそうな値は出たのですが、ただしこのコードまだ未完成です。なぜかというと、センサーアレイはRGGBなのでG(グリーン)が二つあるのですが、その応答が少し違うことに起因します。Gを2つ使って求めると強度の山が2つでできてしまい、ばらつきが大きくなりすぎるからです。そのため今回は2つあるG素子のうち1つだけを使うことにしました。
RGBになりあぷらなーとさんが示してくれた通り、CFAの時よりもばらつきは少なくなることがわかりました。それでもまだノイズが大きすぎるのでCFAの時と同様にunknown factorとしてフィッティングするところに無理やり入れました。RGBの場合、unknown factorはノイズ(標準偏差の２乗)をルート2で割ってやることでZWO、SharpCapの結果とかなり一致することがわかりました。
ASI294MC Proで使ったファイルもアップロードしておきます。よかったら参考にしてください。
ASI294MCPro.zip
今回はまだ苦肉の策で、無理やり答えを合わせるように測定結果を割っているだけなので、考察と呼べるようなものにはなりません。CFAの場合出てきたノイズを半分に、RGBの場合ルート2でわってやるとSharpCapやZWOの結果とかなり近い結果になりましたが、その理由は全然わかっていません。
前回のpythonコードを使った時はCFAやRGBの変換はPixInsightを使いました。でも今回はRAWファイルから直接計算してCFAもRGBも取り出しています。ここら辺はMatlabに移ったことでやりやすくなったところだと思います。このように今回はかなり噛み砕いてブラックボックスがないように進めてきたつもりです。撮影した画像はSharpCapのスクリプトで撮影したものと同等と仮定していいかと思います。これでも結果が一致しないということは、
週末を結構な時間費やしましたが、コンバージョンファクターに関しては結局昨年から進展がなかったので疲れてしまいました。それでも得られた成果としては、
それでもMatlabが使いやすいことは分かったので、今後の画像解析はMatlabで進められそうな目処はついたと思います。
元々この解析はダークノイズ解析につなげたいとか、EOS 6Dのユニティーゲインを測定したいとかで始めたのですが、なかなか進みません。いつ実現できるのかまだまだ未定ですが、諦めることはしたくないのでまたいつか再開する予定です。
私も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でのノイズの少なさには辿りつきません。
いろいろやっていて結局行き着いたところはこの画面で、
小さいと思われる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行目すべてを見たい場合は
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.zip
メーカー値はコンバージョンファクターが3.6程度、ユニティーゲインが110程度ですので、グラフ内の数値と比較してみるとそこそこ一致していることがわかります。
また、以前pythonで書いたコードと比較して同等の結果を出すことできています。
Maltabでもこれまでと同様のことができることがわかり、かつ遥かにシンプルに書けることがわかりました。
カラー版ASI294MC Proの場合: CFAで解析
さて、Matlabでの画像処理にも慣れたところで問題のカラーのASI294MC Proの解析です。
結論から言うと、結局今回も謎は解けずじまいでした。
まずはシンプルなCFAの方から。RAW画像からCFAへ分離するところはビルトインの関数とかはないようなので自分で書いてみました。以前は平均偏差(mad)で無理やり答えを合わせましたが、今回ASI290MMの結果なども含めていろいろ考えた末に、結局諦めて標準偏差(std)を使うことにしました。
各測定点はそのまま解析した結果を表示していますが、フィッティングは無理やり合うように補正してます。考え方としては、謎な部分をunknown factorとしてフィッティングするところで無理やり補正するというやりかたです。今回のCFAの場合、unknown factorはノイズ(標準偏差の２乗)を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つだけを使うことにしました。
RGBになりあぷらなーとさんが示してくれた通り、CFAの時よりもばらつきは少なくなることがわかりました。それでもまだノイズが大きすぎるのでCFAの時と同様にunknown factorとしてフィッティングするところに無理やり入れました。RGBの場合、unknown factorはノイズ(標準偏差の２乗)をルート2で割ってやることでZWO、SharpCapの結果とかなり一致することがわかりました。
考察
ASI294MC Proで使ったファイルもアップロードしておきます。よかったら参考にしてください。
ASI294MCPro.zip
今回はまだ苦肉の策で、無理やり答えを合わせるように測定結果を割っているだけなので、考察と呼べるようなものにはなりません。CFAの場合出てきたノイズを半分に、RGBの場合ルート2でわってやるとSharpCapやZWOの結果とかなり近い結果になりましたが、その理由は全然わかっていません。
前回のpythonコードを使った時はCFAやRGBの変換はPixInsightを使いました。でも今回はRAWファイルから直接計算してCFAもRGBも取り出しています。ここら辺はMatlabに移ったことでやりやすくなったところだと思います。このように今回はかなり噛み砕いてブラックボックスがないように進めてきたつもりです。撮影した画像はSharpCapのスクリプトで撮影したものと同等と仮定していいかと思います。これでも結果が一致しないということは、
- ZWOもSharpCapも何か特殊なことをやっているか
- もしくはまだ私が馬鹿なミスを犯しているか
まとめ
週末を結構な時間費やしましたが、コンバージョンファクターに関しては結局昨年から進展がなかったので疲れてしまいました。それでも得られた成果としては、
- Matlabでの画像解析の手法に慣れることができた。
- あぷらなーとさんが出してくれた、debayerの方がCFAよりもノイズが0.8倍程度と小さくなることを確かめることができた。
それでもMatlabが使いやすいことは分かったので、今後の画像解析はMatlabで進められそうな目処はついたと思います。
元々この解析はダークノイズ解析につなげたいとか、EOS 6Dのユニティーゲインを測定したいとかで始めたのですが、なかなか進みません。いつ実現できるのかまだまだ未定ですが、諦めることはしたくないのでまたいつか再開する予定です。
