目次
CIAOのWavdetect とは
wavdetectは、X線画像の中から局所的に明るい構造を探すためのツールとして知られています。以下の画像のように、実際の研究でも使われている実践的なツールです。
("Probing the nature of the X-ray source IGR J16327-4940 with Chandra"
arXiv:2309.09850より引用。Chandra/ACIS-Iで観測した IGR J16327-4940 周辺のX線画像に、CIAO wavdetect で検出されたX線源候補を重ねた例。)
wavelet は直訳すると「小さな波」で、限られた範囲だけに現れる波のような関数のことです。画像解析ではこの wavelet を画像上で動かしながら、似た形の構造がどこにあるかを調べます。
wavdetect は単に明るいピクセルを探すのではなく、画像の各場所に Mexican Hat wavelet というフィルターを重ねて中心が周囲よりどれくらい盛り上がっているかを調べています。以下はMexican Hat waveletの画像です。
(Cavusoglu (2014), Fig. 1, “Mexican hat wavelet. (a) 3D view. (b) Top view.”より引用。中心が正、周囲が負の形をしており、画像中の局所的な盛り上がりを強調する。)
上の画像で見たように、Mexican Hat wavelet は「中心が正、周囲が負」の形をしています。
この形を画像の各位置に重ねると、中心だけが周囲より明るい場所では正の反応が大きくなり、一様に明るいだけの場所では正負が打ち消されます。
下のGIFでは、その様子を単純化して可視化してみました。GIFはPythonで書いています。
Pythonのコード
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.patches import Circle
from scipy.ndimage import convolve
# ---------------------------------------------------------------------------
# Settings
# ---------------------------------------------------------------------------
np.random.seed(42)
N = 128 # image size [pix]
BG = 2.0 # uniform background [cts/pix]
THRESH = 5.0 # detection threshold [sigma]
SCALES = [1.0, 2.0, 4.0, 8.0]
# Sources: (x_center, y_center, sigma_PSF, total_counts)
# Amplitudes set so each source has SNR ~ 7 near its optimal scale
# sigma_opt = (1 + sqrt(2)) * sigma_PSF ~ 2.414 * sigma_PSF
SOURCES = [
( 32, 32, 1.0, 125.0),
( 32, 96, 2.0, 250.0),
( 96, 32, 4.0, 500.0),
( 96, 96, 8.0, 1000.0),
]
LABELS = ["A", "B", "C", "D"]
plt.rcParams.update({"font.size": 10, "axes.titlesize": 10, "figure.dpi": 110})
# ---------------------------------------------------------------------------
# Core functions
# ---------------------------------------------------------------------------
def make_kernel(sigma, nsig=5):
"""
2D zero-mean Laplacian-of-Gaussian (LoG):
psi(r; sigma) = (1 - r^2 / (2*sigma^2)) * exp(-r^2 / (2*sigma^2))
proportional to -nabla^2 G(r; sigma)
Key properties:
integral(psi) = 0 [zero-mean: local background cancels]
zero-crossing: r = sigma * sqrt(2)
negative minimum: r = 2*sigma, psi_min = -1/e ~ -0.368
integral(psi^2) = pi*sigma^2/2
"""
n = max(5, int(2 * nsig * sigma)) | 1 # odd size
k = np.arange(n) - n // 2
Xk, Yk = np.meshgrid(k, k)
r2 = Xk**2 + Yk**2
s2h = sigma**2 / 2.0 # sigma^2 / 2
return (1.0 - r2 / (2.0 * s2h)) * np.exp(-r2 / (2.0 * s2h))
def wavelet_map(img, sigma):
"""Wavelet coefficient map: W_sigma = f (*) psi_sigma (cross-correlation)"""
return convolve(img.astype(float), make_kernel(sigma), mode="reflect")
def noise_level(sigma, bg):
"""
Standard deviation of wavelet coefficients under Poisson background:
sigma_W = sqrt(bg * integral(psi^2)) = sqrt(bg * pi * sigma^2 / 2)
"""
return np.sqrt(bg * np.pi * sigma**2 / 2.0)
# ---------------------------------------------------------------------------
# Simulate Poisson X-ray image
# ---------------------------------------------------------------------------
x1d = np.arange(N)
XX, YY = np.meshgrid(x1d, x1d)
expected = np.full((N, N), BG, dtype=float)
for (xc, yc, sp, A) in SOURCES:
G = np.exp(-((XX - xc)**2 + (YY - yc)**2) / (2 * sp**2)) / (2 * np.pi * sp**2)
expected += A * G
img = np.random.poisson(expected).astype(float)
# Pre-compute all significance maps
W_maps = []
S_maps = []
sn_vals = []
for sigma in SCALES:
W = wavelet_map(img, sigma)
sn = noise_level(sigma, BG)
W_maps.append(W)
S_maps.append(W / sn)
sn_vals.append(sn)
# Measure peak SNR in aperture around each source (for labels)
APT = 5
peak_snr = np.zeros((len(SCALES), len(SOURCES)))
for i in range(len(SCALES)):
for j, (xc, yc, sp, A) in enumerate(SOURCES):
y0, y1 = max(0, yc - APT), min(N, yc + APT + 1)
x0, x1 = max(0, xc - APT), min(N, xc + APT + 1)
peak_snr[i, j] = S_maps[i][y0:y1, x0:x1].max()
# ---------------------------------------------------------------------------
# Build figure
# ---------------------------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(9.5, 8.5))
fig.subplots_adjust(left=0.07, right=0.97, top=0.91, bottom=0.07,
hspace=0.35, wspace=0.32)
ax_img, ax_ker = axes[0]
ax_W, ax_S = axes[1]
# --- [0,0] Input image (static) ---
im_in = ax_img.imshow(np.log1p(img), cmap="inferno",
origin="lower", extent=[0, N, 0, N])
ax_img.set_title("Input image f(x,y) [log(1+cts)]")
ax_img.set_xlabel("x [pix]"); ax_img.set_ylabel("y [pix]")
fig.colorbar(im_in, ax=ax_img, fraction=0.046, pad=0.04)
for (xc, yc, sp, A), lab in zip(SOURCES, LABELS):
ax_img.plot(xc, yc, "o", color="cyan", ms=5, zorder=5)
ax_img.text(xc + 2, yc + 1, lab, color="cyan", fontsize=9, fontweight="bold")
theta = np.linspace(0, 2 * np.pi, 64)
ax_img.plot(xc + sp * np.cos(theta), yc + sp * np.sin(theta),
color="cyan", lw=0.7, ls="--", alpha=0.6)
# --- [0,1] Kernel (animated) ---
psi0 = make_kernel(SCALES[0])
n0 = psi0.shape[0]
k0 = np.arange(n0) - n0 // 2
ext0 = [k0[0] - .5, k0[-1] + .5, k0[0] - .5, k0[-1] + .5]
im_ker = ax_ker.imshow(psi0, cmap="RdBu_r", origin="lower",
extent=ext0, vmin=-1, vmax=1)
ax_ker.axhline(0, color="gray", lw=0.4, ls="--")
ax_ker.axvline(0, color="gray", lw=0.4, ls="--")
ax_ker.set_xlabel("Dx [pix]"); ax_ker.set_ylabel("Dy [pix]")
fig.colorbar(im_ker, ax=ax_ker, fraction=0.046, pad=0.04)
zero_circle = Circle((0, 0), SCALES[0] * np.sqrt(2),
color="white", fill=False, lw=1.2, ls=":")
ax_ker.add_patch(zero_circle)
title_ker = ax_ker.set_title("")
# --- [1,0] Wavelet coefficients (animated) ---
W0 = W_maps[0]
sn0 = sn_vals[0]
im_W = ax_W.imshow(W0, cmap="RdBu_r", origin="lower",
extent=[0, N, 0, N],
vmin=-4 * sn0, vmax=4 * sn0)
ax_W.set_xlabel("x [pix]"); ax_W.set_ylabel("y [pix]")
cb_W = fig.colorbar(im_W, ax=ax_W, fraction=0.046, pad=0.04)
cb_W.set_label("W_sigma [cts*pix^2]")
title_W = ax_W.set_title("")
cont_W = [None]
snr_texts_W = [ax_W.text(xc + 2, yc + 1, "", color="gray", fontsize=8)
for (xc, yc, sp, A) in SOURCES]
# --- [1,1] Significance map (animated) ---
S0 = S_maps[0]
im_S = ax_S.imshow(S0, cmap="plasma", origin="lower",
extent=[0, N, 0, N], vmin=-2, vmax=10)
ax_S.set_xlabel("x [pix]"); ax_S.set_ylabel("y [pix]")
cb_S = fig.colorbar(im_S, ax=ax_S, fraction=0.046, pad=0.04)
cb_S.set_label(f"S_sigma = W / sigma_W [sigma]")
title_S = ax_S.set_title("")
cont_S = [None]
snr_texts_S = [ax_S.text(xc + 2, yc + 1, "", color="gray", fontsize=8)
for (xc, yc, sp, A) in SOURCES]
main_title = fig.suptitle("", fontsize=13, fontweight="bold")
# ---------------------------------------------------------------------------
# Animation update function
# ---------------------------------------------------------------------------
def update(frame_i):
sigma = SCALES[frame_i]
W = W_maps[frame_i]
S = S_maps[frame_i]
sn = sn_vals[frame_i]
r_zero = sigma * np.sqrt(2) # zero-crossing of LoG
# [0,1] kernel
psi = make_kernel(sigma)
n = psi.shape[0]
k = np.arange(n) - n // 2
ext = [k[0] - .5, k[-1] + .5, k[0] - .5, k[-1] + .5]
im_ker.set_data(psi)
im_ker.set_extent(ext)
ax_ker.set_xlim(k[0] - .5, k[-1] + .5)
ax_ker.set_ylim(k[0] - .5, k[-1] + .5)
zero_circle.set_radius(r_zero)
title_ker.set_text(
f"psi_sigma(x,y) (sigma = {sigma:.0f} pix)\n"
f"zero-crossing r = sigma*sqrt(2) = {r_zero:.1f}, "
f"noise sigma_W = sqrt(bg*pi*sigma^2/2) = {sn:.1f} cts*pix^2"
)
# [1,0] wavelet coefficients
im_W.set_data(W)
im_W.set_clim(-4 * sn, 4 * sn)
if cont_W[0] is not None:
cont_W[0].remove()
if W.max() > THRESH * sn:
cont_W[0] = ax_W.contour(W, levels=[THRESH * sn],
colors="lime", linewidths=1.2,
extent=[0, N, 0, N], origin="lower")
else:
cont_W[0] = None
title_W.set_text(
f"Wavelet coeff. W_sigma = f (*) psi_sigma (sigma = {sigma:.0f} pix)\n"
f"green contour: W_sigma = {THRESH:.0f}*sigma_W = {THRESH*sn:.1f}"
)
for j, (txt, (xc, yc, sp, A), lab) in enumerate(
zip(snr_texts_W, SOURCES, LABELS)):
snr = peak_snr[frame_i, j]
is_det = snr >= THRESH
txt.set_text(f"{lab}\n{snr:+.1f}s")
txt.set_color("lime" if is_det else "gray")
# [1,1] significance map
im_S.set_data(S)
if cont_S[0] is not None:
cont_S[0].remove()
if S.max() > THRESH:
cont_S[0] = ax_S.contour(S, levels=[THRESH], colors="red",
linewidths=1.5,
extent=[0, N, 0, N], origin="lower")
else:
cont_S[0] = None
title_S.set_text(
f"Significance S_sigma = W_sigma / sigma_noise (sigma = {sigma:.0f} pix)\n"
f"red contour: S_sigma = {THRESH:.0f} sigma (detection threshold)"
)
for j, (txt, (xc, yc, sp, A), lab) in enumerate(
zip(snr_texts_S, SOURCES, LABELS)):
snr = peak_snr[frame_i, j]
is_det = snr >= THRESH
txt.set_text(f"{'*' if is_det else ' '}{lab}\n{snr:+.1f}s")
txt.set_color("white" if is_det else "gray")
main_title.set_text(
f"wavedetect — scale sigma = {sigma:.0f} pix "
f"(frame {frame_i+1}/{len(SCALES)})"
)
return [im_ker, im_W, im_S]
ani = animation.FuncAnimation(
fig, update,
frames=len(SCALES), interval=1800,
blit=False, repeat=True,
)
ani.save("wavedetect_2x2_anim.gif",
writer=animation.PillowWriter(fps=0.6), dpi=100)
import matplotlib.pyplot as plt; plt.close("all")
print("Saved: wavedetect_2x2_anim.gif")
左上は仮想の入力画像、右上はそのとき使っている wavelet、左下は画像とwaveletの相関で得られる wavelet coefficient map、右下はそれをノイズで割った significance map です。
左下の wavelet coefficient map は、画像と wavelet の形がどれくらい合っているかを表しています。赤っぽい場所ほど、その scale の Mexican Hat wavelet によく似た「中心が周囲より盛り上がった構造」があることを意味します。一方、青っぽい場所は負の反応で、Mexican Hat の形とは逆向きの構造を表します。
右下の significance map は、wavelet coefficient をノイズで割ったものです。つまり、左下で見えている反応が、ノイズに比べてどれくらい有意かを見ています。この図では、赤い等高線が 5σ の検出しきい値を表しています。
注目したいのは、wavelet scale σ を変えると、強く反応する構造のサイズも変わることです。小さいscaleでは小さなclumpに、大きいscaleではより広がった構造に反応します。つまり wavdetect の scales は、表示の拡大率ではなく、どの大きさの構造を探すかを決めるパラメータです。
詳しく知りたい方はこちら
ここでは、上のGIFで出てきた
- 右上：メキシカンハット wavelet
- 左下：wavelet coefficient map
- 右下：significance map
それぞれが数式で何を表しているのか整理します。
今回はWavdetectについて説明している以下の論文を参考にしました。
"A Wavelet-Based Algorithm for the Spatial Analysis of Poisson Data" (Freeman et al.)
https://arxiv.org/abs/astro-ph/0108429
より正確な内容を知りたい方はご参照ください。ここでは、ざっくりと概要を押さえます。
1. ガウシアンからメキシカンハット waveletを作る
まず、2次元Gaussianを考えます。対称な場合、半径
$$
r^2 = x^2 + y^2
$$
を使って、
$$
G_\sigma(x,y)=\frac{1}{2\pi\sigma^2}
\exp\left(
-\frac{x^2+y^2}{2\sigma^2}
\right)
$$
と書けます。
メキシカンハット waveletは、このGaussianから作られる中心が正、周囲が負のシンプルなwaveletです。Freeman et al. ではメキシカンハット wavelet として、概ね次のような形で出てきます。
$$
\psi_\sigma(x,y)
\propto
\left(
2-\frac{x^2+y^2}{\sigma^2}
\right)
\exp\left(
-\frac{x^2+y^2}{2\sigma^2}
\right)
$$
ここで比例記号 $\propto$ にしているのは、正規化の取り方によって前につく定数が変わるからです。
この式を見ると、中心 $r=0$ では正、ある半径より外側では負になっていて十分遠くではgaussianが0に落としてくれる、という形になっていることがわかります。
実際、
$$
2-\frac{r^2}{\sigma^2}=0
$$
となる位置は、
$$
r=\sqrt{2}\sigma
$$
です。
つまり、$r<\sqrt{2}\sigma$ では中心の正の領域、$r>\sqrt{2}\sigma$ では周囲の負の領域になります。GIF右上の中心が赤っぽく、周囲が青っぽい構造に対応しています。
2. なぜ一様背景が消えやすいのか
メキシカンハット waveletの重要な性質は、全体を足し合わせるとほぼ0になることです。
$$
\iint \psi_\sigma(x,y),dx,dy = 0
$$
これは、中心の正の部分と周囲の負の部分が打ち消し合うという意味です。
そのため、画像が一様背景だけなら、
$$
D(x,y)=B
$$
として、
$$
\iint \psi_\sigma(x,y)B,dx,dy=
B\iint \psi_\sigma(x,y),dx,dy=
0
$$
となります。
つまり、メキシカンハット waveletは一様に明るい場所にはあまり反応せず、中心が周囲より盛り上がった場所に反応しやすいということになります。
3. wavelet coefficient mapとは何か
GIFの左下に出ている wavelet coefficient map は、入力画像とwaveletの相関を取ったものです。
2次元画像を $D_{i,j}$、waveletを $W$ とすると、Freeman et al. の式では次のように書かれます。
$$
C_{i,j}=\sum_{i'}\sum_{j'}
W_{i-i',j-j'}D_{i',j'}
$$
これは、やっていることは重み付きの和です。
各ピクセル $(i,j)$ を中心としてwaveletを置き、その周囲の画像カウント $D_{i',j'}$ にwaveletの重み $W_{i-i',j-j'}$ を掛けて、全部足しています。
$D_{i,j}$ はそのピクセルに入ったphoton countですが、$C_{i,j}$ はそのピクセル周辺の構造が、scale $\sigma$ のMexican Hat waveletにどれくらい似ているかを表す量です。
4. GIFの左下は何を表示しているのか
上のGIFでは、左下に
$$
W_\sigma = f \ast \psi_\sigma
$$
のような量を表示しています。
ここで、
- $f(x,y)$：入力画像
- $\psi_\sigma(x,y)$：scale $\sigma$ のメキシカンハット wavelet
- $W_\sigma(x,y)$：wavelet coefficient map
です。
つまり、GIF左下は入力画像をscale $\sigma$ のメキシカンハット waveletで見たときの反応マップということになります。
赤っぽい場所は、画像とwaveletの形がよくあっているということを表しています。
つまり、そのscaleで見ると中心が周りよりも盛り上がった構造があることを表現しています。
青っぽい場所は負の反応です。 これは、メキシカンハットの形とは逆向きの構造、つまり中心よりも周りの方が相対的に強いような構造に対応しています。
実際に検出してみる
実際の天体画像でWavedetectを試してみましょう。天体のデータはChandra Data Archiveから探します。
Target Nameで天体の名前から検索することもできます。今回はCassiopeia Aと検索して、このデータで解析してみたいと思います。
カシオペアAについては、以下の記事に詳しいです。
CIAOというソフトと、ds9というソフトを使います。
インストールについては公式情報を参照してください。
一連の流れについて、酒井さんのこちらの記事が参考になります。
ChandraアーカイブのObsID 9117には、すぐに使える全体画像（Level 2プロダクト）が primary/ ディレクトリに入っています。
${DS9} ${DATA_DIR}/primary/acisf09117N003_full_img2.fits.gz &
DS9が起動したら、メニューから初期設定をします。
| 設定項目 | 操作 |
|---|---|
| カラーマップ | Color → Heat |
| スケール | Scale → Log |
| ズーム | Zoom → Fit |
DS9のFile → Open のGUIファイル選択は .fits.gz を認識できず、ファイルが全て灰色になります。ターミナルからの起動かFinderからのドラッグ＆ドロップで開いてください。
画像上で右クリックしたままマウスを左右にドラッグすると、輝度のMin/Maxをリアルタイムで変更できます。
次に、作業ディレクトリの作成とCIAOの初期化を行います。
mkdir -p ${DATA_DIR}/Mywork
cd ${DATA_DIR}/Mywork
source ${CIAO_DIR}/bin/ciao.sh
CIAO 4.18.0 ... と表示されれば初期化完了です。
wavdetectの実行には、カウント画像と露光マップの2つが必要です。fluximage がイベントファイルから両方を一括で生成してくれます。
fluximage \
infile="${DATA_DIR}/primary/acisf09117N003_evt2.fits.gz[EVENTS]" \
outroot="casa9117" \
bands="broad" \
asol="${DATA_DIR}/primary/pcadf09117_000N001_asol1.fits.gz" \
badpixfile="${DATA_DIR}/primary/acisf09117_000N003_bpix1.fits.gz" \
binsize=4 \
psfecf=0.393 \
clobber=yes
| パラメータ | 値 | 説明 |
|---|---|---|
bands |
broad |
0.5–7.0 keV（Chandra CSC定義のbroad band） |
binsize |
4 |
1ピクセル≈1.97″ |
psfecf |
0.393 |
PSFエネルギー集中率のデフォルト値 |
実行には少し時間がかかります。完了後、以下のファイルが生成されます。
casa9117_broad_thresh.img ← wavdetectの入力画像
casa9117_broad_thresh.expmap ← 露光マップ
casa9117_broad_thresh.psfmap ← PSFマップ
casa9117_broad_flux.img ← フラックス画像
ファイルが揃っているか確認します。
ls casa9117_broad_thresh.img casa9117_broad_thresh.expmap
ここから、実際にwavdetectで点源検出を実行します。
wavdetect \
infile="casa9117_broad_thresh.img" \
outfile="casa9117_src.fits" \
scellfile="casa9117_scell.fits" \
imagefile="casa9117_imgfile.fits" \
defnbkgfile="casa9117_nbkg.fits" \
expfile="casa9117_broad_thresh.expmap" \
psffile="casa9117_broad_thresh.psfmap" \
scales="1.0 2.0 4.0 8.0 16.0" \
sigthresh=1e-6 \
clobber=yes
| パラメータ | 値 | 説明 |
|---|---|---|
scales |
"1.0 2.0 4.0 8.0 16.0" |
探索する構造の空間スケール。小さいほど細かい構造を検出できる |
sigthresh |
1e-6 |
検出の有意性しきい値 |
参考 :
最後に、ds9に検出結果を重ねます。
${DS9} \
casa9117_broad_thresh.img \
-scale log \
-cmap heat \
-regions load casa9117_src.fits &
実行結果
2007年のカシオペアAで実行してみると、次のような画像を得ます。
上の画像は bands="broad"（0.5–7.0 keV）、scales="1.0 2.0 4.0 8.0 16.0"、binsize=4（1ピクセル≈2″）での wavdetect 結果です。
緑の楕円が検出された領域で、大小さまざまなサイズが混在しているのがわかります。
Cas Aは点源ではなく広がった天体なので、SNR内の構造が片っ端から検出されています。
これは、wavdetectが周囲より局所的に明るい領域を拾った結果です。
このままでは大きな拡散構造と細かい点状構造が混在しているため、次のステップでは scales を絞り込んで細かい構造だけを見ていきます。
上の画像は scales="1.0 2.0" に絞り込んだ結果です。前の全スケール結果と見比べると変化が明確です。
| 比較項目 | 全スケール（1.0〜16.0） | 細かいスケール（1.0〜2.0） |
|---|---|---|
| 楕円のサイズ | 大小混在 | 小さく均一 |
| 内部の拡散構造 | 大きい楕円で覆われる | ほぼ消える |
| 外縁のフィラメント | 大きい楕円が連なる | 小さい楕円が点在 |
| 検出の性質 | 拡散構造 ＋ 点状構造 | 比較的コンパクトな局所構造が中心 |
細かい構造だけ見たい場合は小さい scales に絞る、SNR全体の広がりを把握したい場合は大きい scales も加える、という使い分けがあります。
参考
緑の楕円は、wavdetect が検出した候補です。Cas A内部には knot やフィラメントが多いため、これらをすべて本物の点源とみなすことはできません。一方で、今回の結果では、Holland-Ashford et al. (2024) が位置合わせに用いた校正星も検出できています。同論文では、wavdetect で検出した候補を Gaia DR3 と照合し、複数 epoch での一致や SNR 本体の放射から分離できることを確認した上で、校正用点源として採用しています。
(Holland-Ashford et al. (2024) より引用。)
そのため、この結果はCas A内部の局所構造も拾っているが、先行研究で使われた信頼できる点源も含めて検出できていると解釈できます。
実行結果 : SS 433の場合
CasAでは、広がった超新星残骸内部のフィラメントなども wavdetect に反応することを見ました。次に、明るい点源を含む別の観測例として SS433 を見て、検出器由来の構造が wavdetect の結果にどう影響するかを観察してみたいと思います。
1. 生データをDS9で確認する
イベントファイルをそのまま DS9 で開くとこのようになっています。
中央の白輝点が SS 433。縦に伸びる線は後述する readout streakというもの。中心のSS 433が非常に明るく、周囲に背景天体と思われる点源が散らばっているのがわかります。
2. wavdetect をかけると検出器の構造が見えてしまう
FOV フィルタをかけてイメージを作成し、wavdetect を実行します。
dmcopy "acisf03790N004_evt2.fits.gz[EVENTS][sky=region(acisf03790_000N004_fov1.fits.gz)]" evt_fov.fits
dmcopy "evt_fov.fits[bin x=::4,y=::4]" img_bin4_fov.fits
wavdetect infile=img_bin4_fov.fits outfile=src_fov.fits scellfile=scell_fov.fits \
imagefile=imgout_fov.fits defnbkgfile=nbkg_fov.fits \
scales="1 2 4 8 16" sigthresh=1e-6 psffile=none clobber=yes
緑のリージョンが視野の四辺に集中しています。
この結果を見ると、検出領域が視野端やチップ境界付近にも多く現れています。これは、必ずしもそこに物理的な天体があるという意味ではありません。ACIS-S画像ではReadout streak などが画像上に強い構造として現れることがあります。Readout streak についての詳しい内容はこちらが参考になります。
3. フィルタ処理を整えてもう一度実行する
以下の 3 点を改善して再実行します。
- エネルギーフィルタ：0.5–8 keV に絞る
- ビンサイズ：bin=4 → bin=2
- Readout streak のマスク：SS 433 の明るさで生じる縦の線を除外
※ 除外前の画像。中央の線はReadout streakである可能性が高いです。
dmcopy "acisf03790N004_evt2.fits.gz[EVENTS][sky=region(acisf03790_000N004_fov1.fits.gz)][energy=500:8000]" evt_filtered.fits
dmcopy 'evt_filtered.fits[sky=region(streak_only.reg)][bin x=::2,y=::2]' img_bin2_masked.fits clobber=yes
wavdetect infile=img_bin2_masked.fits outfile=src_masked.fits scellfile=scell_masked.fits \
imagefile=imgout_masked.fits defnbkgfile=nbkg_masked.fits \
scales="1 1.414 2 2.828 4" sigthresh=1e-6 psffile=none clobber=yes
streak_only.reg は DS9 で描いた box リージョンで、field() で全域を選択し -box(...) でストリーク部分を除外しています。
今回はreadout streakの影響を見るために簡略化してpsffile=noneとなっていますが、厳密な点源のカタログを作る場合はPSFマップの設定が必要です。
※ DS9で描いたboxリージョン
physical
field()
-box(4034.9509,4086.6092,57.60001,1020.0002,356.00283)
最終的に以下のような画像を得た。
チップ端の偽検出が排除され、視野内の点源（背景 AGN・X線星など）が適切に検出されました。最終的な検出数は68個です。
4.HRCの観測結果を見る
比較のために、Chandra/HRCで取得されたSS433の画像も見てみます。
HRCはACISとは異なり、CCDではないため、ACISで問題になったような readout streak は基本的に現れません。そのため、明るい中心源から縦方向に伸びる人工的な線状構造はかなり抑えられます。
一方で、HRC画像では検出器の有効視野、露光端、背景イベントの分布が画像上に強く見えることがあります。したがって、HRCで wavdetect を実行した場合も、検出された領域をそのまま物理的な点源とみなすのではなく、検出器の視野形状や背景構造の影響を考慮する必要があります。
ACISでは readout streak やチップの構造、HRCでは検出器の有効視野や背景分布が、それぞれ wavdetect の結果に影響します。wavdetect は「天体かどうか」を判断しているのではなく、画像上で メキシカンハット wavelet に強く反応する source-like structure を拾っている、という点が重要です。
Wavdetectの原理
Wavdetectの原理は次の論文が参考になります。
wavdetect の一番中心にある計算は、次の式です。
$$
C_{i,j} = \sum_{i'}\sum_{j'}
W_{i-i',j-j'}D_{i',j'}
\tag{7}
$$
$$
C_{i,j}
\equiv
\langle W \star D\rangle_{i,j}
\tag{8}
$$
難しく見えますが、やっていることは 重み付き足し算 です。
ここで、
- $\ D_{i,j}$：入力画像のピクセル(i',j') のカウント数
- $\ W_{i-i',j-j'}$：中心を ((i,j)) に置いたwaveletの重み
- $\ C_{i,j}$：その位置でのwavelet coefficient
です。
つまり、各ピクセル (i,j) について、その場所を中心にMexican Hat waveletを置き、周囲の画像カウントと掛け算して足し合わせるということをしています。
もっと単純化して、1次元で考えてみます。
たとえば、Mexican Hatのミニ版として
$$
W = [-1,\ 2,\ -1]
$$
という重みを考えます。
これは中心が正、周囲が負の形というMexican Hat Waveletを再現しています。
ここで一様背景を導入します。
$$
D = [5,\ 5,\ 5]
$$
として、
$$
C = (-1)\times5 + 2\times5 + (-1)\times5 = 0
$$
となります。
一方、中心だけ明るければ、
$$
D = [5,\ 12,\ 5]
$$
なので、
$$
C = (-1)\times5 + 2\times12 + (-1)\times5 = 14
$$
となり、大きな正の値が残ります。
ここで注意したいのは、棒グラフの高さは入力カウント$\ D$ を表している、ということです。一方、棒の上に書いた値は、各位置のカウントに wavelet の重みを掛けた寄与 $\ W\times D$ を表しています。つまり、中心のカウントが 12 で、そこに wavelet の重み $+2\ $が掛かるので、
$$
12 \times (+2) = +24
$$
となります。
左右のカウントはそれぞれ 5 で、wavelet の重みは (-1) なので、
$$
5 \times (-1) = -5
$$
です。
最後にこれらの寄与をすべて足すと、
$$
C = (-1)\times 5 + 2\times 12 + (-1)\times 5 = 14
$$
となります。
つまり、wavdetect の第一段階では、画像全体にこのようなフィルタを走査して、どこがMexican Hat的な構造に見えるかを表す wavelet coefficient map を作っています。
この一連を動画にしてみました。
これは、カーネルをスライドして$\ C_i $ を作っているところです。最初のGIFで見た wavelet coefficient map がどのように作られているかが見えてくると思います。入力画像をコピーしているわけではなく、あくまで局所的な盛り上がりを見ています。
ここで得られる$\ C_{i,j}$ は、元画像のカウント数そのものではありません。
$$
C_{i,j} \neq D_{i,j}
$$
$\ D_{i,j}$ は、そのピクセルに入ったphoton countです。 一方、$\ C_{i,j}$ はそのピクセル周辺の構造が、指定したscaleのMexican Hat waveletにどれくらい似ているかを表します。
そのため、全体的に明るい領域であっても、一様に明るいだけなら $\ C_{i,j}$ は大きくなりません。逆に、背景の中に小さな盛り上がりがあれば、そこでは$ \ C_{i,j}\ $が大きくなります。
この意味で、wavdetect は「明るい場所を探す」というより、周囲に対して中心が盛り上がった source-like structure を探すアルゴリズムだと言えます。
フーリエ変換との関連性
wavelet変換は、フーリエ変換を知っている人には「使う波を変えた変換」と考えるとわかりやすいです。
フーリエ変換では、信号を
$$
e^{ikx}
$$
のような、空間全体に広がった波に分解します。
そのため、「どの周波数成分がどれくらい含まれているか」を調べるのが得意です。
一方、wavelet変換では、マザーウェーブレットと呼ばれる局所的な波を用意し、それを平行移動・拡大縮小しながら信号や画像に重ねます。
つまり、フーリエ変換が信号を「周波数ごとの波」に射影する変換だとすれば、wavelet変換は信号を「位置とスケールを持った小さな波」に射影する変換です。
wavdetect では、このマザーウェーブレットとして Mexican Hat wavelet を使っています。画像の各位置に Mexican Hat wavelet を重ねて、さらに scale を変えながら相関を計算します。これによって、「どこに、どの空間スケールの source-like structure があるか」を探すことができます。
左のフーリエ基底は「どの周波数が含まれるか」を見るのに向いていて、右の Mexican Hat wavelet は「どの位置に、どのスケールの局所構造があるか」を見るのに向いています。フーリエ変換では周波数$\ k$を変えるのに対して、wavelet変換では位置とscaleを変えます。
フーリエ変換は「全体に含まれる周波数」を見るのが得意で、wavelet変換は「局所的な構造がどこにあるか」を見るのが得意。 wavdetect の scales は、どの大きさの構造を見るかに対応している。
wavdetect は単に明るいピクセルを探すコマンドではありません。メキシカンハットwavelet を画像上で動かして、scale を変えながら周囲に対して局所的に明るい構造を探すアルゴリズムです。
そのため、点源だけでなくCas Aのような広がった天体のフィラメント等にも反応します。一方で、readout streak や検出器由来の構造にも反応するので、検出結果をそのまま天体とみなすのではなく、画像や検出器構造と照らし合わせて解釈することが重要です。