DoGフィルタによる星検出

星検出とは

星検出とは画像内の明るい点を特定し、それらを星として検出するプロセスです。天体写真などのスタッキンング処理などで、星の位置を特定するために使用されます。この記事ではDoGフィルタを用いた星検出アルゴリズムについて解説します。

DoGフィルタ

DoG(Difference of Gaussian)は2つのわずかに異なる標準偏差を持つガウシアンフィルタを用いて、画像のスケール空間での点状構造を検出する手法。LoG(Laplacian of Gaussian)と似ていますが、計算コストが低く、スケール空間での点状構造の検出に適しています。簡単にいうとぼかしの違いを取ることで「小さな明るい点」を強調する方法で、星検出においては、星のような小さな明るい点を背景から際立たせるのに有効です。

ガウス関数

ガウス関数（Gaussian function）とは、最も基本となる正規分布の形状を持つ曲線を表します。2次元のガウス関数は以下の式で表されます。

ガウスフィルタ

ガウスフィルタ(Gaussian filter)はガウシアンぼかしとも呼ばれ、画像の各ピクセルを周辺ピクセルの重み付き平均に置き換える処理のことです。ガウス関数に従って計算された重み（カーネル）を周囲のピクセルに適用し、中心からの距離が近いピクセルの値をより強く反映させます。これにより、画像のノイズを減らし、滑らかなぼかし効果を得ることなどができます。ただし、ガウスフィルタはエッジをぼかすため、シャープな輪郭を保ちたい場合、メディアンフィルタなどの他のフィルタと組み合わせることが推奨されます。

DoGの星検出におけるアルゴリズム

画像に対して2種類のぼかしをかけます

• 小さいぼかし（σ_small）
• 大きいぼかし（σ_large）

そして、それらの差を取ります

ここで、星と星以外の反応の違いは以下のようになります。

• 星（小さい明るい点）: DoGが大きい
  • 小ぼかしでは残る
  • 大ぼかしでは消える
• 星以外（広い構造）: DoGが小さい
  • 両方で似たようにぼける

DoGはLoGの近似となっています。

RustでのDoGフィルタ実装例

このとき、以下のようになります。

• x：中心からの距離
• σ：広がり（ぼけの強さ）
• exp(⋅)：自然指数関数（eのべき乗）

ガウシアンカーネル

カーネル関数とは、データの類似度を測るための関数で、ここではガウシアンカーネルを用いて画像の各ピクセルを周囲のピクセルの重み付き平均に置き換えるための重み（カーネル）として使用します。 ガウシアンカーネルのことをRBF kernelとも呼ぶ。データ同士の距離が近いほど値が大きくなり類似度が高いことを示します。

カーネル半径

ガウス関数

正規化係数

ゆえに正規化後のガウシアンカーネルは以下のようになります。

ここで、正規化後のガウシアンカーネルは、以下を満たします。

fn gaussian_kernel_1d(sigma: f32) -> Vec<f32> {
    let radius = (sigma * 3.0).ceil().max(1.0) as i32;
    let mut kernel = Vec::with_capacity((radius * 2 + 1) as usize);
    let mut norm = 0.0_f32;
 
    for i in -radius..=radius {
        let x = i as f32;
        // ガウス関数
        let v = (-x * x / (2.0 * sigma * sigma)).exp();
        kernel.push(v);
        norm += v;
    }
 
    for k in &mut kernel {
        *k /= norm;
    }
 
    kernel
}

ガウスフィルタの実装

2次元ガウシアンは $G(x,y)=Gx(x)Gy(y)$という性質があります。そのため、1次元ガウシアンカーネルを縦と横に1次元畳み込みで計算することで計算量を $O(N^2)$ から $O(2N)$ に削減することができます。 $\sigma = 0$の時は何もしません。

fn gaussian_blur(img: &GrayImageF32, sigma: f32) -> Vec<f32> {
    if sigma <= 0.0 {
        return img.data.clone();
    }
 
    let kernel = gaussian_kernel_1d(sigma);
    let radius = (kernel.len() / 2) as i32;
 
    let mut temp = vec![0.0_f32; img.data.len()];
    let mut out = vec![0.0_f32; img.data.len()];
 
    for y in 0..img.height {
        for x in 0..img.width {
            let mut sum = 0.0_f32;
            for (k_idx, k) in kernel.iter().enumerate() {
                let dx = k_idx as i32 - radius;
                sum += img.get_clamped(x as i32 + dx, y as i32) * k;
            }
            temp[(y * img.width + x) as usize] = sum;
        }
    }
 
    for y in 0..img.height {
        for x in 0..img.width {
            let mut sum = 0.0_f32;
            for (k_idx, k) in kernel.iter().enumerate() {
                let dy = k_idx as i32 - radius;
                let sy = (y as i32 + dy).clamp(0, img.height as i32 - 1) as u32;
                sum += temp[(sy * img.width + x) as usize] * k;
            }
            out[(y * img.width + x) as usize] = sum;
        }
    }
 
    out
}

ここでget_clamped関数は、画像の境界を超えた場合に端のピクセル値を返すようにしています。

fn get_clamped(&self, x: i32, y: i32) -> f32 {
    let cx = x.clamp(0, self.width as i32 - 1) as u32;
    let cy = y.clamp(0, self.height as i32 - 1) as u32;
    self.data[self.index(cx, cy)]
}

DoGによる星検出

以下を満たす画素を星候補として検出する。

1. DoGが閾値 $\theta$ 以上である (だいたい $0.01 \lt \theta \lt 0.10$ )
2. SNRが $\tau$ 以上である
3. 線状エッジ上の反応を除去
4. 8近傍平均より十分に突出している
5. 8近傍で局所最大である

DoGが閾値 $\theta$ 以上である

閾値を上げると検出数は減るが、誤検出は減りやすいです。 だいたい $0.01 \lt \theta \lt 0.10$ に設定するのが良いです。

SNRが$\tau$以上である

SNR(Signal-to-Noise Ratio)とは、信号とノイズの比率を示す指標であり、局所標準偏差で正規化した値。

線状エッジ上の反応を除去

線状エッジ上の反応を除去し、角・点状ピークを優先します。線状エッジ上の反応は、DoGの値が大きくても、星ではない可能性が高いためです。

まず、DoGの2次微分(Hessian)を近似します。 DoGを $D(x,y)$ とすると、中心 $(x,y)$ で次を計算します。

これで Hessian 行列を作ります。

ここで$\det(H)\le 0$（実装では微小値以下）の場合は鞍点や不安定な形状で、点状ピークらしくないため、除外します。

エッジ判定スコアRを以下のように定義します。

固有値を $\lambda_1,\lambda_2$ とすると

固有値比を $r=|\lambda_1/\lambda_2| \ (r\ge1)$ とおくと

ここで、Rの値は以下のように解釈できます。

• $R$ が大きい: 一方向だけ曲率が強い・$r$ が大きい・
• $R$ が小さい: 2方向の曲率が近い・$r \approx 1$・点状ピーク

閾値比を $r_{\max}$ として以下を満たせば、点状ピークとして通過させます。

ゆえに、Rの値が大きい場合は線状エッジの可能性が高く、除外されやすいです。逆に、Rの値が小さい場合は点状ピークの可能性が高く、通過しやすいです。

8近傍平均より十分に突出している

8近傍平均との差（prominence）が小さい弱いピークを除外します。

画素座標を $(x,y)$、DoG画像を $D(\cdot,\cdot)$ とします。

まず、8近傍集合を以下のように定義する。

ここで、prominenceは以下の式で表せます。

• $P(x,y)>0$ : 中心が周囲平均より高い（山になっている）
• $P(x,y)\approx 0$ : 周囲とほぼ同じ高さ
• $P(x,y)<0$ : 中心が周囲平均より低い（谷になっている）

したがって、突出係数を $\alpha$、閾値を $\varepsilon$ として、以下の条件を満たす場合にのみ、星候補として採用します。

8近傍で局所最大である

画素座標を $(x,y)$、DoG画像を $D(\cdot,\cdot)$ とします。

まず、8近傍集合を以下のように定義します。

以下の条件を満たす場合にのみ、中心が厳密な局所最大となり、星候補として採用します。

誤検知の除外

NMSの除外

NMS(Non-Maximum Suppression)とは非極大抑制のことで、局所的に最大値でない候補を除外する手法です。ここでは、星候補の中で、候補同士が近すぎる場合は応答の強い方を残すようにします。これにより、同じ星が複数の候補として検出されることを防ぎます。

候補点を $c_i = (x_i, y_i, r_i)$とし、最小距離を$d_{\min}$とします。採用集合を $S$ とすると、$c_{\pi(k)}$ を見て、次を満たすときだけ採用します。

すなわち、どの採用済み点からも距離が $d_{\min}$ より大きい場合のみSに追加し、それ以外は除外します。

孤立点の除外

孤立点を誤検知とみなし除外します。周囲に他の星候補が存在しない場合は、ノイズや誤検知の可能性が高いです。

2点 $d_i, d_j$ のユークリッド距離二乗は以下のようになります。

$d_j$ が $d_i$ の近傍である条件は、半径を $R$とすると、

ゆえに、各点 $d_i$ の近傍数は以下のように書けます。

最低必要近傍数を $K$とすると、以下を満たすもの以外は除外できます。

実際の検出結果

実際に以下の画像を用いて、DoGフィルタを用いた星検出を行います。

出力結果は以下のようになりました。黄色の+マークが検出された星の位置を示しています。

参考文献

1. 【画像処理】空間フィルタリング_ラプラシアン/LoG/DoGフィルタ - Deliberate Learning
2. 画像処理の基礎 ― OpenCV による基本的な例 ―
3. ガウス過程回帰〜カーネルについてわかりやすく説明〜