生画像へのマスクの掛け方と測光する画像についての議論

[COIAS-program]

事象

COIASサーバの /diskCOIAS/data04/ecliptic/9960/9960-3,4/2019-11-01 に置いてある5枚の画像について解析を行ったところ、マスクなし画像では移動天体が写っているのに、マスクあり画像では消えているものがあることが分かった。 以下の添付スクショはその一例で、3枚目のx=1432,y=740付近で確認できる。 上のスクショはマスクなし、下のスクショはマスクありのもので、どちらも3枚目の同じ場所を撮ったものである。 赤四角で囲まれた部分について、マスクなし画像では光源があるのに、マスクあり画像では光源が消えている。

この光源は左上から右下に移動している光源であり、3枚目以外には写っていないので、マスクがかかるべきではない。 ただし、このピクセル座標の値には1枚目はNaNが入っており、5枚目は画像が欠けていて負の値が入っている。

[COIAS-program]

直接的原因

直接的な原因はマスクの掛け方のアルゴリズムと画像の悪さにある。マスクは web-coias-back-app/src2-startsearch2R/subm2.py で使用されているアルゴリズムと同じもので掛けている。

まずsubm2.pyの段階で解析するfits画像の構造とマスクの掛け方を説明する。 hdul = fits.open([fits画像名]) でfits画像を開いてデータをhdulに格納すると、hdul[0].dataには生データが、hdul[1].dataにはマスク用データが入る。hdul[0].dataもhdul[1].dataも2次元配列で、hdul[0].data[y][x]に(x,y)のピクセルの値が入っている。生データは文字通り生データで、明るければ大きい値が入るし、負の値が入ったりNaNが入ったりもする。 一方でマスク用データは、16bitのunsigned intであり、暗ければ小さい値が、明るければ大きい値が入り、基本的に暗い場合は0が入る。 subm2.pyでは、5枚の画像のマスク用データを使用し、あるピクセル座標において得られる5つのマスク用データの中央値が0でない場合には、そのピクセルにマスクをかける、という処理を行う。

以下の添付スクショは、上記「事象」で問題を確認したx=1432,y=740における5枚の画像の生データの値とマスク用データの値を調べたものである(y座標が2100-740になっているのは、png座標とfits座標で上下反転しているから。2100はy方向の総ピクセル数)。 この画像を見ると確かにこの座標における5つのマスク用データのうち3つがノンゼロであり、中央値もノンゼロなのでマスクをかける条件に抵触していることがわかる。 しかし2個おかしい点があって、

• 1枚目の生データにNaNが入っている場合のマスク用データの値が256である
• 5枚目の生データに負の値が入っている場合のマスク用データの値がゼロではない

NaNピクセルに大きいマスク用データが入っているのは、NaNのピクセルにマスクをかけて欲しいという意図のためだと思われる。

[COIAS-program]

間接的影響

マスクをかけるアルゴリズムは配布されているものではなく、自前で書いたものであるため問題が起きる可能性がある。 さらに上記のように変なデータが入っている時には、その可能性はさらに高くなる。 現状では、ある位置にある天体の明るさ(フォトン数=天体が写っている範囲内の生データのカウント数の合計)を測定する際に、マスクありデータを用いている。 マスクの掛け方が信用ならない以上、それを用いて測光を行うことは危険かもしれない。

そのため、測光を行う際はマスクなし画像を使った方が良いかもしれない。

[COIAS-program]

直接的原因の解決案

直接的原因で書いたように、生データがNaNのピクセルに256というマスク用データが入っていることが問題の原因の1つになっている。 また生データが負の値なのにマスク用データにゼロではない値が入っているのも原因の1つだろう。 ただし、上記2点はこのピクセルでそうだということを確認したのみで、一般的であるかどうかは分からない。

上記を鑑みて、以下の2つの解決案が思いつく。この2つは排他的なものではなく、両方使ってもいいかもしれない。

1. 生データがNaNのピクセルのマスク用データには強制的に0を入れてしまう。
2. マスクを掛けるか否かの中央値の値に閾値を設ける。現状では中央値が0でない場合にはマスクをかけるとしているが、例えば中央値が6より大きい場合にはマスクをかける、などの閾値を設ける。もしそうすれば、直接的原因で書いた事例では中央値が5であるので、このピクセルにはマスクがかからないことになる。

[COIAS-program]

直接的原因の解決案の懸念点

しかし上記の直接的原因の解決案を何も考えずに実行すれば良いかというと疑問が残る。本来マスクをかけるべき場所にマスクがかからないようになったり、別の影響が出る恐れがある。 1.の解決案の場合、例えば5枚の画像のあるピクセルにおいて、2枚に恒星が写っており、3枚がNaNの場合、このピクセルにはマスクがかからないことになるが、これは望ましい結果ではない。現状のアルゴリズムではこのピクセルには正しくマスクがかかる。 2.の解決案の場合、そもそも閾値の値をいくつにすれば良いか議論しなくてはいけない。

いずれにしても、マスクを掛けるアルゴリズムを修正すると、全ての画像に影響してしまう。そのため様々な画像を用いて、上記の解決案をオリジナルCOIASにて色々試してみて、慎重に検討した方が良い。

[COIAS-program]

間接的影響の解決案

マスクなし画像に対して測光を行えば良さそうに思える。

ただし、マスクなし画像に対して光源検出(SExtractor)および検出した光源を用いた移動天体候補自動検出を行ってしまうと、恒星ばかり検出してしまって移動天体の自動検出に大きな影響を及ぼしてしまう。 そのため、光源検出と自動検出 <= マスクあり画像、測光 <= マスクなし画像、の使い分けが必要である。

現状で解析を行なっている画像は、ビニング後にマスクをかけ、マスクとNaNの領域をスカイのノイズで埋めたものを使用している(PytyonがNaNを踏むとエラーで止まる可能性があるため)。 しかしマスクをかけていないが、NaNはスカイのノイズで埋めた画像は用意していないので、別途それを用意する必要がある。

2023/4/8現在、COIAS_program_githubリポジトリのsugiura-nonmask-photometryブランチでは、上記2点の変更が反映されており、マスクなし画像に対して測光を行うようになっている。

[COIAS-program]

間接的影響の解決案の懸念点

マスクなし画像を用いると、移動天体の周囲にいる恒星が消されないため、バックグラウンドを測定する時に恒星の明かりを巻き込んでしまいうまく測光できない場合がある。

移動天体の明るさを測るときはその天体の範囲内のフォトン数の合計を求めれば良いが、その時にバックグラウンドも巻き込んでしまうため、バックグラウンドのフォトン数を引き去るという操作を行う。 そのバックグラウンドのフォトン数は、移動天体の範囲の少し外側の領域のフォトン数の平均から見積もる。

移動天体の周囲に恒星があった場合、バックグラウンドの見積もりに恒星を巻き込んでしまうためうまく測定できない。 一方でマスクあり画像の場合は、恒星はマスクで消されているので問題はない。

実例を挙げる。添付のスクショはオリジナルCOIASであるマスクなし画像をGUIで手動測定した時のスクショである。 赤い四角で囲まれている光源が移動天体で、その左上に写っている光源が恒星である。 この移動天体について、sugiura-nonmask-photometryブランチのコードでマスクなし画像に対して測光を行ったところ、恒星のせいでバックグラウンドを大きく見積りすぎてしまい、移動天体のフォトン数 - バックグランドのフォトン数がマイナスとなり、測光できなかった。 一方でmainブランチのコードでマスクあり画像に対して測光をすると、正しく測光ができる。

マスクなし画像は確かに生データに近いので安心かもしれないが、落とし穴は存在するようである。

[COIAS-program]

両解決案の最大の懸念点

現状ではCOIASサーバに置いてある画像は事前にビニングマスクをして、ビニングマスク済み画像も同じ場所に置くことにしているが、そこのアルゴリズムを変えるなら全ての画像に対してビニングマスクをやり直さないといけない。それにはおそらく10日ほどかかる。 さらに、現状ではHDDの容量の8割程度まで画像を詰め込んでいるが、アルゴリズムを変更して新たなビニングマスク済み画像を保存する必要になったりするとおそらく溢れる。上書きするなどして使用容量を大きく増やさないようにしないといけない。

5月中にリリースする予定ならあと50日ほどしかないので、何回もアルゴリズムを変えて全てビニングマスクをやり直している時間はない。妥当なアルゴリズムを2週間程度で見つけて、1回だけ全ての画像にビニングマスクをかけ直す時間しかないと思われる。 リリース後にアルゴリズム更新をしても良いが、その場合はビニングマスクを完了させるために10日ほどサービスを停止することになる。

別の考え方として、事前のビニングマスクを諦めるという手もある。 この場合、サーバにビニングマスク後の画像を置く必要がなくなるので、使用容量が半分くらいまで減るというメリットがある。 一方でユーザが解析する時にビニングマスクを実行しなければいけなくなるので、探索準備モードにおける解析時間が少なくとも30秒以上は伸びてしまい、ユーザビリティに悪影響を及ぼす。

[kn1cht]

@COIAS-program 追いつきました．分かりやすくまとめてくださり，ありがとうございます．

1. hdul[1].data に入るマスク用データは使いにくい点があるということですが，これは我々ではなくFITS画像を配布する側が用意したものですか？　それとも， subm2.py より前のどこかで生成されるものですか？
2. マスク用データが使いにくいということであれば，それを使わずに生データを直接加工して（NaNや負の値をゼロに置換するなど）マスク計算に使うのではだめでしょうか？
3. 直接的原因の解決案の懸念点1について：これは突き詰めていくと，移動天体を見逃すか，恒星を誤検出するかのトレードオフになりそうに感じます（たとえば，2枚でもピクセルが光っていればマスクをかけるという基準に変えれば2枚恒星，3枚NaNでもマスクがかかるが，大きめの小惑星が通過中で連続して同じ位置が光っているという場合でもマスクが適用されてしまう）．最後の方で仰ったように，解析画像を選んだときにどれだけNaNが混入してしまうのかという点は要確認だと思います．

[COIAS-program]

@kn1cht ご確認いただきありがとうございます。

1. hdul[1].data に入るマスク用データは使いにくい点があるということですが，これは我々ではなくFITS画像を配布する側が用意したものですか？　それとも， subm2.py より前のどこかで生成されるものですか？

subm2.pyにてhdul[1].dataに入っている配列はFITS画像配布者が用意したものです (あまり詳しくないですが、すばる望遠鏡の生画像を一次処理した時に生成されるものかと。いずれにしても、アーカイブに上がっている画像の時点でhdul[1].dataは用意されています)。正確には、アーカイブデータの段階ではhdul[2].dataにマスク用データが格納されていますが、binning.pyで2x2にビニングしたのちにhdul[1].dataとして格納しています。そういう意味では元画像のマスク用データにビニングするという処理はしていますが、元画像のhdul[2].dataも32bit unsigned intなので、ビニング以外の処理はしていません。

2. マスク用データが使いにくいということであれば，それを使わずに生データを直接加工して（NaNや負の値をゼロに置換するなど）マスク計算に使うのではだめでしょうか？

ここは僕からはなんとも。。。 @urakawa-seitaro いかがでしょうか?

3. 直接的原因の解決案の懸念点1について：これは突き詰めていくと，移動天体を見逃すか，恒星を誤検出するかのトレードオフになりそうに感じます（たとえば，2枚でもピクセルが光っていればマスクをかけるという基準に変えれば2枚恒星，3枚NaNでもマスクがかかるが，大きめの小惑星が通過中で連続して同じ位置が光っているという場合でもマスクが適用されてしまう）．最後の方で仰ったように，解析画像を選んだときにどれだけNaNが混入してしまうのかという点は要確認だと思います．

全くおっしゃる通りで、結局トレードオフになると思います。なので色々試して最適だと思われる条件 (閾値とかNaNの扱いとかマスクの基準とか) を決めていくしかないですね。。。

[kn1cht]

@urakawa-seitaro @COIAS-program あれからいろいろと調べていましたら，HSCデータ解析講習という資料を見つけました． https://hsc.mtk.nao.ac.jp/pipedoc/DataWorkshop2022/_downloads/2682edf0a494c8b005acf4efa51952ae/hsc_kaiseki_20221128.pdf

p. 92に，HDU2に入ったMask Imageについての説明があります． この説明を見る限り，マスク用データに入っている値は十進の数字に意味があるわけではなく，ビットマスクなのではないでしょうか？　例えば，5が入っている場合，0101となるので1番目と3番目のフラグが立っている…という感じです．

また，上でNaNが256になる話がありましたが，256は2進数で 100000000 なので，このスライドでPlane8（9番目）がNO_DATAを意味しているのと合います．

[kn1cht]

なお，ヘッダのマスクレイヤーを参照と書いてあるので手元のサンプルwarp画像をfv等で開いてみましたが，該当しそうな表記は見つけきれませんでした

[kn1cht]

The Hyper Suprime-Cam Software Pipeline (Publ. Astron. Soc. Japan (2014) 00(0), 1–39) のTable 2に，BAD, SAT, INTRP, CR, ...の意味の説明がありました． 総合すると，何かしらよくないピクセルであるときに0以外の値が入るということなので，天体が写っているかどうかの判定をマスク用データが非ゼロかどうかで行うのがそもそも間違いかもしれません

[COIAS-program]

@kn1cht @urakawa-seitaro 上の伊東さんの解釈は合っていそうですね。 注意としては、資料に「ビットとマスクイベントの対応は不定」とあるように、資料に乗っているどのフラグがどの桁に対応しているのかは変わりうる点でしょうか。 実際に適当なfitsを開いて、hdul[2].headerを見てみたら、フラグは13種類あり、9にBRGHT_OBJECT、12にUNMASKED_NANが入っていました。

他のトレードオフ要素と違いマスク用データの解釈違いは致命的なので、以下のようにbinning.pyとsubm2.pyを書き換えた上で、COIASサーバにある全画像のビニングマスクをやり直すことを提案します。(もちろん事前にテストはしますが) この時、ついでにマスクなし画像で測光できるように、使えるマスクなし画像も生成しておきます。

• binning.pyにてhdul[2].headerを確認して、DETECTEDとDETECTED_NEGATIVEフラグが立つ2進数での桁を確認して書き出す
• binning.pyでビニングする時、マスク用データ(hdul[2].data)は平均を取らず、左下のピクセルの値を採用する。(nbin=2の場合はビニングする2x2の4ピクセルのうち左下、nbin=4なら4x4の一番左下。どこを選ぶかという自由度はあるが適当で良いと思われる。)
• subm2.pyでマスクを掛けるアルゴリズムを以下のように変更:
  • DETECTEDかDETECTED_NEGATIVEフラグが立っているピクセルには1を、そうでないピクセルには0をセットする
  • 測定画像枚数でそのmedianを取り、それがノンゼロならマスクを掛ける。

これでマスク用データの本質を捉えたマスク処理ができるはずです。 上記「直接的原因の解決案の懸念点」での問題は解決できないですが、あまり気にしなくても良いのかなと思います。

上記「間接的影響」に関連して1つ思い出したことがあります。 現状では4枚の画像の組が多く、しかもそれなりな割合で極端に画像間隔が狭いものがあったりします。 そうすると移動天体であっても星像の1部が2枚の画像で被ったりするので、その場合マスク処理の条件に抵触して星像が欠けたり最悪消えたりします。(だいたいの場合は欠けるくらいで済むと思いますが) 欠けるだけだったらマスクあり画像でも目で見て移動天体があると分かるので手動測定できますが、マスクあり画像に測光すると測光値がおかしくなってしまいます。 なので、浦川さんが主張されるようにマスクなし画像への測光が安全かなと思います。 この時「間接的影響の解決案の懸念点」が問題となってくるので、昨日の伊東さんの提案通りバックグラウンドの値にスカイの平均を使うようにする。

ひとまず上記全ての変更をオリジナルCOIASに加えてみて、天体が消えた画像を含めいろいろな画像でテストを行い、変更前と後でどれだけ良くなったか・見た目に問題はないか・測光した値は大きく変わらないか、などをチェックしてみようと思います。

[kn1cht]

ありがとうございます．対策について同意です．

[COIAS-program]

参考: HSC pipeline manual https://hsc.mtk.nao.ac.jp/pipedoc/pipedoc_4/j_tool/index.html

[COIAS-program]

＊sugiura-nonmask-pohotometryのskyの明るさは移動天体周りのスカイで推定。 ＊移動天体探しは、新しく修正したマスク画像で行う。ただし測光はマスクなし画像に対して行う。 ＊マニュアルモードで、四角形を作る時に表示されているのはマスク画像。ただし測光と位置測定はマスクなし画像に対して実行。

上記にて議論はまとまり、オリジナルCOIAS・ウェブCOIASの両バックエンド側リポジトリに反映が完了したので、closeします。