先日のエントリで11月30日早朝に撮った写真にボリソフ彗星(2I/Borisov)が写っているのか写っていないのか結論が出せなくてもやもやしていました。
なんか心霊写真みたいにそう見てしまうとそう見えてしまうという現象という可能性が捨てきれず、何らかの客観的な指標を見て判断したいのだが… というところに落ち着いたわけですが、あの後色々解析めいたことをやってみました。
方法
p値の計算
彗星と思われる光点(以下光点)が背景光のノイズではない(可能性が高い)ということを示すため、光点の周辺 320 x 320 ピクセルの領域(背景領域)のピクセル輝度の平均値が、光点の周辺 3 x 3 ピクセルの領域の輝度の平均値以上である確率(p値)を片側t検定で求めました。p値が十分小さければノイズではないと考えられます。
厳密には背景領域からは恒星の写っている部分を取り除くべきですが、適切な方法を思いつかなかったためそのまま平均値や分散を計算しています。これにより平均値も分散も大きくなるため、p値は実際より大きくなる(よりノイズと判定されやすくなる)はずです。
p値の視覚化
比較のため光点の場所だけでなく、画像周囲の1ピクセル幅の領域を除く全ての点についても同様にp値を計算し、p値が特定の値以下の領域を視覚化した画像を作成しました。
p < 0.001, p < 0.0001, p < 0.00001, p < 0.000005, p < 0.000003 を視覚化した画像を半透明化して重ね合わせることで5段階のグレースケール(p 値が小さいほど白い)によるヒートマップも作成しました。
光点の位置の確認
p値のヒートマップに Stellarium で表示した星図を重ね合わせることで、光点に対応する p 値のピークの位置と星図上の彗星の位置が一致するかどうかを目視で確認しました。
星図の位置合わせのため、恒星基準コンポジットした画像をレベル補正して星を強調し、星図のサイズに合わせて301%拡大したものを用意しました。
恒星の位置がなるべく一致するように Stellarium で表示した星図を重ねたものが以下です。
星図の 2I/Borisov の 04:48 (コンポジット画像の reference frame の撮影時刻)の位置はヒートマップを重ねた時にも視認できるように赤くしてあります。また彗星周囲の光芒はヒートマップと重なると輝度が違って見えるため消してあります。
対象となる画像
処理対象となる画像は彗星核基準コンポジット後の画像ですが、画像全体を処理するとデータ量が膨大になるため、また画像の一部に見られる顕著な傾斜カブリの影響を避けるため、光点を中心とした 320 x 320 ピクセルの領域を切り抜いて使用しました。
DeepSkyStacker でスタックしたTIFF画像は 32bit/チャンネル の形式ですが、画像処理で使用するライブラリ(rmagick)ではピクセル値を 16bit 整数値しか取れなかったので Photoshop CC で切り抜きした画像を保存する際に 16bit/チャンネル の形式に変換しました。
処理プログラム
ruby による処理プログラムを作成しました。画像の入出力には rmagic を、t検定の計算には statsample を使用しました。ソースコードは付録に掲載します。
結果
p値が特定の値以下の領域は以下のようになりました。表示は2倍に拡大しています。
ノイズである確率が0.1%以下の領域です。面積は4344ピクセル(4.24%)。恒星の光跡がはっきり見えます。画像左上には傾斜カブリによるノイズらしきものが見られます。中央の光点周辺は面積21ピクセルの独立した領域になっています。
ノイズである確率が0.01%以下の領域です。面積は1478ピクセル(1.44%)。光点周辺は面積11ピクセルの独立した領域になっています。
ノイズである確率が0.001%以下の領域です。面積は401ピクセル(0.39%)。恒星の光跡がだいぶかすれてきました。光点周辺は面積4ピクセルの独立した領域になっています。
ノイズである確率が0.0005%以下の領域です。面積は252ピクセル(0.25%)。恒星の光跡がだいぶかすれてきました。光点周辺は面積3ピクセルの独立した領域になっています。
ノイズである確率が0.0003%以下の領域です。面積は170ピクセル(0.17%)。恒星の光跡がだいぶかすれてきました。光点周辺の面積は1ピクセルだけになりました。この点のp値は 0.00000285 です。約35万分の1の確率です。
これらの画像を重ね合わせてヒートマップにしたものを恒星と星図の画像に重ねたものが以下になります。
2I/Borisov の位置に光点がほぼ重なることがわかります。
以下は3倍に拡大したものです。
よく見ると光点のp値の最小値の位置(光点の中心位置)と星図の 2I/Borisov の位置との間にはズレがあります。光点の中心位置は 2I/Borisov から5ピクセル左、4ピクセル上にずれています。上の画像は3.01倍したものを3倍しているので、等倍の画像上では0.55ピクセル左、0.44ピクセル上にズレていることになります。
考察
光点がノイズである確率は十分低いように思われます… 最初は p < 0.01 で有意とか思っていたのですが、それだとどう見てもノイズにしか見えない部分が大量に出てきたのでした。
よく考えると10万ピクセルあるので1000ピクセルくらいはノイズが浮いてくるはずです。が、それよりもずっと多い10922ピクセル(10.7%)が該当してしまいました。p < 0.001 でもかなりノイズが見えるし、有意水準をどのくらいにするのが妥当なのかよくわからなくなってきました。
光点の中心がノイズの確率が0.00000285ということは同じレベルのノイズが350877ピクセルに1つ出るということです。このレベルのノイズが 2I/Borisov の位置を含む350877ピクセルの中に発生した場合、 2I/Borisov の位置と一致する確率を考えます。
光点の中心は 2I/Borisov の位置から左に0.55ピクセル、上に0.44ピクセル離れた場所にあり、その距離は0.70ピクセルです。半径0.70ピクセルの円の面積は1.54ピクセルです。350877ピクセルの領域にランダムに発生したノイズが 2I/Borisov の位置を中心とする半径 0.70 ピクセルの領域に入る確率は 1.54 / 350877 = 0.00000439 です。
こんな計算でいいんでしょうか？もっとも、可視化した結果を見ると実際にノイズが発生する確率はp値より二桁ぐらい高いようにも見えますが… とはいえ二桁違ってもまだ0.05%以下です。
また、前回のエントリで「これが彗星だと言うのなら、この写真彗星だらけでは？」と言いましたが、p < 0.000003 の画像を見ると光点の他はほぼ恒星の光跡部分しか残っていないように見えます。
光点が何らかの天体であるとして 2I/Borisov の位置からわずかにズレているのは気になります。星図との重ね合わせの精度が気になりますが恒星の位置ズレを目視で数点確認した限りでは等倍で0.5ピクセル以内のようです。光点の位置ずれはそれより40%以上大きいので何か理由がありはずです。
ズレが発生する要因としては彗星核基準コンポジットの際の彗星の位置指定のズレが考えられます。実際の彗星の軌道のベクトルが指定した始点と終点を結ぶベクトルと一致しなければ、彗星像はブレた形でスタックされ、輝度のピークはズレてしまいます。
DeepSkyStacker の Info ファイルに指定した彗星位置の座標が記録されているので、そこからベクトルがわかります。それによると始点は左上原点のピクセル座標系で (2512.00 , 1924.00) 終点は (2496.00, 1965.00) でした。
しかし、ガイドエラーで先頭フレームと最終フレーム(= reference frame)の写野はズレていて、先頭フレームの dX は 1.92、dY は 1.12 でした。dX, dY はフレーム上の星像をそれだけ動かすと reference frame 上の星像に重なるという量です。なので、先頭フレームに付けた彗星マークの位置 (x, y) を最終フレームの座標系に直すと (x + dX, y + dY) になり、(2512.00 + 1.92, 1924.00 + 1.12) = (2513.92, 1925.12) になります。
よって、先頭フレーム上の彗星は (2496.00 - 2513.92, 1965 - 1925.12) = (-17.92, 39.88) だけ移動して最終フレームにスタックされます。つまり、最終フレーム上の彗星の位置から右に 0.02, 下に 0.68 ピクセルするズレてスタックされます。彗星像が対称に広がっていると仮定するとスタックされた像のピーク位置はその中間点になるので、右に 0.01、下に 0.38 ピクセルずれるはずです。
って、あれ？縦方向のズレが実際と逆ですね？実際に写った光点のピーク位置は、本来写るはずの位置から (-0.55 + 0.01, -0.44 + 0.38) = (-0.56, -0.78) だけズレていることになります。ズレが大きくなってしまいました… ズレの距離は0.96ピクセルになりました。
p値のヒートマップに本来写るはずの彗星像の位置を示したのが以下の図です。水色の点で写るはずの彗星像を模式的に示しています。
改めてノイズがこのズレの範囲内に出る確率を計算します。半径0.96ピクセルの円の面積は2.90ピクセル。2.90 / 350877 = 0.00000827 です。
というわけで位置ズレについては謎が残ってしまいましたが、確率的にノイズでないのはほぼ確実で、このあたりに光跡を残すような十分明るい恒星も見当たらないのでこれがボリソフ彗星だと思うことにします…
そもそも統計とか専門でないのでこんな計算でいいのかもよくわからないんですが、そうだと思うことにします！
付録: 処理プログラムのソースコード
ライブラリとして rmagic と statsample が必要です。どちらも ruby gem として入手可能です。
p値の計算
p値の計算は重いので、表をテキストデータとして出力しておいて、視覚化の際に利用します。
#!/usr/bin/ruby # coding: utf-8 # USAGE: p-map.rb image-file # 周辺 3x3 ピクセルの p 値の表を出力する require 'statsample' require 'rmagick' abort "USAGE: p-map.rb image-file" if ARGV.length != 1 img = Magick::Image.read(ARGV[0]).first width = img.columns height = img.rows sky = Array.new(width * height) sky_i = 0 sky_map = Array.new(height) { Array.new(width) } img.each_pixel { |px, x, y| pi = px.intensity() sky_map[y][x] = pi sky[sky_i] = pi sky_i += 1 } sky_v = sky.to_vector sky_mean = sky_v.mean sky_sd = sky_v.sd sky_n = sky_v.size star = Array.new(3 * 3) for y in 0..(height - 1) do for x in 0..(width - 1) do if y < 1 || height - 2 < y || x < 1 || width - 2 < x then print "#{1.0}\t" next end i = 0 for sy in (y-1)..(y+1) do for sx in (x-1)..(x+1) do star[i] = sky_map[sy][sx] i += 1 end end star_v = star.to_vector star_mean = star_v.mean star_sd = star_v.sd star_n = star_v.size t = Statsample::Test::T.two_sample_independent(star_mean, sky_mean, star_sd, sky_sd, star_n, sky_n, false) df = Statsample::Test::T.df_not_equal_variance(star_sd, sky_sd, star_n, sky_n) p = Statsample::Test.p_using_cdf(Distribution::T.cdf(t, df), :right) print "#{p}\t" end puts end
p値の視覚化
p-map.rb で出力した表からp値が指定した値以下の領域を示す画像を出力します。
#!/usr/bin/ruby # coding: utf-8 # USAGE: map-to-image.rb in-map out-image threashold require 'rmagick' include Magick abort "USAGE: map-to-image.rb in-map out-image threashold" if ARGV.length != 3 map_file = ARGV[0] image_file = ARGV[1] t = ARGV[2].to_f map = [] f = File.open(map_file) f.each { |line| row = line.chomp.split("\t") values = Array.new(row.size) row.each_with_index { |s,i| values[i] = s.to_f } map.push(values) } w = map[0].size h = map.size img = Image.new(w, h) map.each_with_index { |values, y| values.each_with_index { |v, x| color = v < t ? "white" : "black" img.pixel_color(x, y, color) } } img.write("#{image_file}")
