FITS 画像の smoothing - 宇宙線実験の覚え書き

複数の望遠鏡の角度分解能、beam size の違いを吸収するためや、少数統計の画像を見やすくする目的で、FITS 画像に smoothing をかけることがたまにあります。自分が使っている手段と、その注意点を簡単にまとます。

DS9 を使った smoothing

次のように command line から DS9 を引数付きで起動すれば、kernel radius が 20 pixel のgaussian で smoothing してくれます。

$ ds9 image.fits -smooth yes -smooth radius 20

ここで、"kernel radius" の定義に注意が必要です。20 pixel の gaussian で smoothing をかけると言った場合、通常は σ = 20 pixel だと考えます。しかし DS9 の場合、σ = (kernel radius)/2 が定義です*1。したがって、上記の command で得られる画像は、σ = 10 pixel でぼかしたものになります。Kernel radius には整数値しか使えないので、用途が限定的になります。

同等の機能は、GUI の menu から Analysis > Smooth Parameters... と辿ることで使うこともできます。また、ぼかしたものを FITS に保存することも可能です。

XIMAGE を使った smoothing

HEAsoft が入っていれば、XIMAGE が一緒に入っているはずです。DS9 でやったのと同じ作業をしたければ、次のようにして下さい。FITS が整数値の画像である場合、0 未満の値が切り捨てられるので、そういう場合には scaling_factor という引数が必要です。

[XIMAGE> read/fits image.fits
[XIMAGE> smooth/sigma=20/scaling_factor=1000
[XIMAGE> disp
[XIMAGE> write/fits smooth.fits

DS9 とは異なり、sigma の値に実数を使うことが可能です。

注意点

恐らく両者とも内部で FFT の演算をしています。そのためか、gaussian の裾の計算は、ある程度の範囲で省略してしまいます。DS9 では kernel radius よりも外の pixel については、ぼかしの処理を行いません (つまり、半径 2σ の範囲でしか計算を行わない)。また、XIMAGE では最大で 50 × 50 pixel の範囲内でしか計算しません (これは σ の値に依り、小さめの σ を与えると 30 × 30 pixel だったりします)。

裾を切り捨てた場合、積分値はどうなるのか気になるところですが、積分値は正しい値が保持されます。つまり、裾を切り捨てた分だけ、peak が少し高くなります。

DS9 と XIMAGE で kernel radius、もしくは sigma の値を 20 に設定した場合の出力結果を比較したものが、以下の図です。元画像は、201 × 201 pixel の、中央の pixel のみが 255、他が 0 の値を持つ FITS 画像です。DS9 のほうが、広がりが半分しかないことが分かります。また、DS9 では裾が半径 20 pixel で切り捨てられ、XIMAGE では 50 × 50 pixel で切り捨てられているのが分かります。黒の lego plot が FITS 画像、赤線が領域を限って fitting した gaussian です。

参考までに、作図に使った python script は以下の通りです*2。

import ROOT
import fitsroot

ROOT.gStyle.SetOptFit()
ROOT.gStyle.SetOptStat()

can = ROOT.ExactSizeCanvas("can", "can", 1000, 500)
can.Divide(2, 1, 1e-10, 1e-10)

img = [fitsroot.FitsImage("ds9.fits", 0, "DS9"),
       fitsroot.FitsImage("ximage.fits", 0, "XIMAGE")]

f = [ROOT.TF2("f0", "[0]*exp(-((x - 101)**2. + (y - 101)**2.)/(2.*[1]**2))", 101 - 10, 101 + 10, 101 - 10, 101 + 10),
     ROOT.TF2("f1", "[0]*exp(-((x - 101)**2. + (y - 101)**2.)/(2.*[1]**2))", 101 - 10, 101 + 10, 101 - 10, 101 + 10)]

for i in range(2):
    can.cd(i + 1)
    f[i].SetParameter(0, 10)
    f[i].SetParameter(1, 10)
    img[i].Fit(f[i], "", "goff")
    img[i].GetXaxis().SetRangeUser(101 - 40, 101 + 40)
    img[i].GetYaxis().SetRangeUser(101 - 40, 101 + 40)
    img[i].Draw("surf3")

    maximum = img[i].GetMaximum()*1.1
    minimum = img[i].GetMinimum()
    img[i].SetMaximum(maximum)
    f[i].SetMaximum(maximum)
    
    f[i].SetLineColor(2)
    f[i].SetRange(101 - 40 - 0.5, 101 - 40 - 0.5, minimum, 101 + 40 + 0.5, 101 + 40 + 0.5, maximum)
    f[i].Draw("same surf")