photutils

Photutils - 天体測光

Photutils は、天体画像から定量データを抽出するための現代的な標準です。SExtractor や DAOPHOT などの旧来のツールを置き換え、Astropy エコシステムと統合された、モジュール式で Python らしいインターフェースを提供します。

使用する場面

• 天体画像内の恒星や銀河を検出する場合（DAOStarFinder、IRAFStarFinder）。
• 円形、楕円形、長方形の開口を使って天体の明るさを測定する場合。
• PSF フィッティング測光を行う場合（密集した恒星領域で不可欠です）。
• 複雑な 2D 背景（空の背景）を推定して差し引く場合。
• 天体ソースの重心（正確な中心）を計算する場合。
• 「等輝度解析」（銀河の形状測定）を行う場合。
• しきい値処理に基づいて画像をセグメンテーションする場合。

参考ドキュメント

公式ドキュメント：https://photutils.readthedocs.io/
GitHub：https://github.com/astropy/photutils
検索パターン：photutils.aperture、photutils.detection、photutils.background、photutils.psf

基本原則

まず背景を処理する

背景を正しく差し引かなければ、正確な測光はできません。Photutils は、全体および局所的な空のレベルを推定するツールを提供します。

開口測光と PSF 測光

• 開口測光：固定形状の範囲内にあるピクセルを合計します。孤立した天体ソースや銀河に最適です。
• PSF 測光：恒星の数学的モデルをデータにフィッティングします。重なり合う恒星や高精度な測定に最適です。

単位と WCS

Photutils は単位を完全に認識します。画像に WCS メタデータが含まれている場合、開口はピクセルだけでなく、天球座標（度）で定義できます。

クイックリファレンス

インストール

pip install photutils

標準インポート

import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u
from photutils.detection import DAOStarFinder
from photutils.aperture import CircularAperture, aperture_photometry
from photutils.background import Background2D, MedianBackground

基本パターン - 検出と開口測光

import numpy as np
from photutils.detection import DAOStarFinder
from photutils.aperture import CircularAperture, aperture_photometry

# 1. Detect sources
# fwhm: Full Width at Half Maximum of the stars
finder = DAOStarFinder(fwhm=3.0, threshold=5.*std)
sources = finder(image - background)

# 2. Define apertures at detected positions
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))
apertures = CircularAperture(positions, r=4.0)

# 3. Perform photometry
phot_table = aperture_photometry(image - background, apertures)
print(phot_table)

重要なルール

✅ 推奨事項

• まず背景を差し引く - 背景を差し引いた画像では、検出と測光の精度が常に向上します。
• Background2D を使用する - 空のレベルに変動がある画像では、単一の中央値よりも 2D メッシュによる推定のほうが正確です。
• ゲインを指定する - 誤差（不確実性）を正確に推定するには、aperture_photometry に検出器のゲインを渡します。
• FWHM を合わせる - 検出器の fwhm パラメーターが、画像の実際の PSF と一致していることを確認します。
• 開口をベクトル化する - 1 つの CircularAperture オブジェクトに数千個の位置を格納できます。ループ内でオブジェクトを作成しないでください。
• SkyAperture を使用する - WCS がある場合は、赤経／赤緯で開口を定義し、画像の回転や再投影後も対象を捉えられるようにします。

❌ 非推奨事項

• マスクされたピクセルを無視しない - 画像に不良ピクセル（宇宙線、飽和した恒星）がある場合は、測光関数にマスクを渡します。
• 小さすぎる開口を使わない - 開口が小さすぎると光を取りこぼし（開口補正が必要になります）、大きすぎるとノイズが増えます。
• デフォルトを無条件に信頼しない - 開口を画像に重ねて必ず可視化し、恒星に合っていることを確認します。
• ピクセル単位と天球座標単位を混在させない - (x, y) と (RA, Dec) のどちらを渡すかに注意します。

アンチパターン（絶対に避けること）

# ❌ BAD: Manual pixel summation
# light_sum = np.sum(image[y-r:y+r, x-r:x+r]) # No circular masking!

# ✅ GOOD: Proper aperture photometry
from photutils.aperture import CircularAperture
ap = CircularAperture((x, y), r=r)
phot = aperture_photometry(image, ap)

# ❌ BAD: Detecting sources on raw data with high background
# finder = DAOStarFinder(threshold=100)(raw_image)

# ✅ GOOD: Detection on noise-subtracted data
from photutils.background import Background2D, MedianBackground
bkg = Background2D(image, (50, 50), filter_size=(3, 3), bkg_estimator=MedianBackground())
sources = finder(image - bkg.background)

背景推定（photutils.background）

堅牢な 2D 背景マッピング

from photutils.background import Background2D, MedianBackground
from astropy.stats import SigmaClip

sigma_clip = SigmaClip(sigma=3.0)
bkg_estimator = MedianBackground()

# box_size: size of the mesh grid
bkg = Background2D(image, box_size=(50, 50), 
                   filter_size=(3, 3), 
                   sigma_clip=sigma_clip, 
                   bkg_estimator=bkg_estimator)

print(f"Global median background: {bkg.background_median}")
clean_image = image - bkg.background

ソース検出（photutils.detection）

恒星検出アルゴリズム

from photutils.detection import DAOStarFinder, IRAFStarFinder

# DAOStarFinder: Optimized for stars
daofind = DAOStarFinder(fwhm=3.0, threshold=5.0*std)
sources = daofind(data)

# IRAFStarFinder: Alternative implementation
iraffind = IRAFStarFinder(threshold=5.0*std, fwhm=3.0)
sources_iraf = iraffind(data)

# Extract coordinates for photometry
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))

開口測光（photutils.aperture）

複数の形状の扱い

from photutils.aperture import CircularAperture, CircularAnnulus, aperture_photometry

# 1. Main aperture and Annulus (for local background)
positions = [(100.5, 120.3), (250.1, 88.4)]
aperture = CircularAperture(positions, r=5.0)
annulus = CircularAnnulus(positions, r_in=10.0, r_out=15.0)

# 2. Local background subtraction
# Calculate median in the annulus
annulus_masks = annulus.to_mask(method='center')
bkg_median = []
for mask in annulus_masks:
    annulus_data = mask.multiply(image)
    annulus_data_1d = annulus_data[mask.data > 0]
    bkg_median.append(np.median(annulus_data_1d))

# 3. Perform photometry and subtract local background
phot = aperture_photometry(image, aperture)
phot['local_bkg'] = bkg_median
phot['aper_sum_corrected'] = phot['aperture_sum'] - (bkg_median * aperture.area)

PSF 測光（photutils.psf）

重なり合うソースのフィッティング

from photutils.psf import BasicPSFFitter, IntegratedGaussianPRF
from photutils.background import MMMBackground

# 1. Define the PSF model (Gaussian or custom)
psf_model = IntegratedGaussianPRF(sigma=1.5)

# 2. Initialize the fitter
fitter = BasicPSFFitter(psf_model, fitshape=(11, 11), 
                        finder=daofind, 
                        bkg_estimator=MMMBackground())

# 3. Perform fitting
# result_table contains fitted x, y, and flux
result_table = fitter(image)

# 4. Get the residual image (original - fitted stars)
residual_image = fitter.get_residual_image()

重心の特定（photutils.centroids）

恒星中心の精密な特定

from photutils.centroids import centroid_com, centroid_2dg

# Center of mass (simple)
x1, y1 = centroid_com(image[y-5:y+6, x-5:x+6])

# 2D Gaussian fit (most precise for stars)
x2, y2 = centroid_2dg(image[y-5:y+6, x-5:x+6])

実践的なワークフロー

1. 画像全体のパイプライン（検出から等級まで）

def extract_magnitudes(image, zero_point=25.0):
    # 1. Subtract 2D Background
    bkg = Background2D(image, box_size=(64, 64))
    data_sub = image - bkg.background
    
    # 2. Detect sources
    finder = DAOStarFinder(fwhm=2.5, threshold=5*bkg.background_rms_median)
    sources = finder(data_sub)
    
    # 3. Aperture Photometry
    pos = np.transpose((sources['xcentroid'], sources['ycentroid']))
    ap = CircularAperture(pos, r=4.0)
    phot = aperture_photometry(data_sub, ap)
    
    # 4. Convert to instrumental magnitudes
    phot['mag'] = -2.5 * np.log10(phot['aperture_sum']) + zero_point
    return phot

2. 残差による測光品質の確認

def check_psf_quality(image, sources):
    """Fits PSF and returns residual to check for artifacts."""
    # ... setup PSF fitter ...
    result = fitter(image, init_guesses=sources)
    residual = fitter.get_residual_image()
    
    # If residual has 'holes' or 'donuts', PSF model is bad.
    return residual

パフォーマンスの最適化

フィッティングでの n_jobs の使用

現在の Photutils は、すべてのフィッターで n_jobs をネイティブにサポートしているわけではありませんが、画像をタイルに分割することで並列処理できます。

開口マスク

非常に大きな画像では、画像全体に対して計算する代わりに、aperture.to_mask() を使って関連する部分配列だけを抽出します。

よくある落とし穴と解決策

log10 の負のフラックス

背景差し引きが過剰な場合、暗い星の一部で aperture_sum が負になることがあります。

# ✅ Solution: Filter results before log
valid = phot[phot['aperture_sum'] > 0]
valid['mag'] = -2.5 * np.log10(valid['aperture_sum'])

混雑した視野

星が重なっている場合、開口測光は機能しません。

# ❌ Problem: One aperture contains 3 stars
# ✅ Solution: Switch to PSF Photometry (photutils.psf)
# This fits the shapes and accounts for overlapping light.

座標の順序

photutils では (x, y) を使用し、これは (列, 行) に対応します。

# ❌ Error: image[x, y] # This is wrong for numpy!
# ✅ Solution: Use image[y, x] or image[row, col]

Photutils は、天体画像を測定するための「物差し」です。堅牢な統計的検出と精密な数学的フィッティングを組み合わせることで、天文学者はデジタル画像から、不確実性が明らかな科学的カタログを作成できます。