치수 계산하기 - robust edge detection 2

Phase-Based Edge Detection

이미지의 edge는 일반적으로 위상(phase)의 급격한 변화에서 발생합니다. 특히 Fourier나 Hilbert Transform을 통해 구한 위상 정보를 분석하면 에지 위치를 subpixel 정밀도로 파악할 수 있습니다.

Fourier Transform 기반 방식
- 이미지를 주파수 영역으로 변환 (FFT)
- 특정 방향의 주파수 성분을 필터링
- 위상 정보 추출
- 위상 변화가 급격한 부분을 에지로 판단

Hilbert Transform 기반 방식
- Hilbert Transform은 1D 시그널에서 즉각적인 위상(instantaneous phase) 추정에 사용됩니다.
- 이미지를 방향성 필터와 함께 처리하여 각 방향의 위상 응답을 얻고 에지를 검출합니다.

 

Hilbert Transform은 무엇인가?
Hilbert Transform은 1D 실수 신호를 복소수 analytic signal로 확장해주는 변환입니다. 이 변환을 통해:
- 위상(phase): 실수/허수 쌍으로부터 추출 가능
- Envelope (진폭): 에너지 강도 해석 가능

C++에는 Hilbert Transform이 표준으로 제공되지 않기 때문에 FFT 를 기반으로 구현이 필요하다.

 

 

예제

#include <iostream>
#include <vector>
#include <cmath>
#include <fftw3.h>

std::vector<double> hilbert_phase(const std::vector<double>& signal) {
    int N = signal.size();
    std::vector<double> phase(N);

    // FFTW 메모리 할당
    fftw_complex *in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

    // 신호 복사
    for (int i = 0; i < N; ++i) {
        in[i][0] = signal[i]; // real
        in[i][1] = 0.0;       // imag
    }

    // Forward FFT
    fftw_plan p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(p);

    // Hilbert Transform filter in frequency domain
    for (int i = 1; i < N / 2; ++i) {
        out[i][0] *= 2;
        out[i][1] *= 2;
    }
    for (int i = N / 2 + 1; i < N; ++i) {
        out[i][0] = 0;
        out[i][1] = 0;
    }

    // Inverse FFT
    fftw_plan q = fftw_plan_dft_1d(N, out, in, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_execute(q);

    // 위상 계산
    for (int i = 0; i < N; ++i) {
        double real = in[i][0] / N;
        double imag = in[i][1] / N;
        phase[i] = atan2(imag, real);
    }

    // 자원 해제
    fftw_destroy_plan(p);
    fftw_destroy_plan(q);
    fftw_free(in);
    fftw_free(out);

    return phase;
}

int main() {
    // 간단한 신호 예시
    std::vector<double> signal = {0, 1, 2, 5, 10, 5, 2, 1, 0};
    auto phase = hilbert_phase(signal);

    std::cout << "Phase:" << std::endl;
    for (double p : phase)
        std::cout << p << " ";
    std::cout << std::endl;

    return 0;
}