OpenCV로 배우는 영상 처리 및 응용/ 변환영역 처리

  • 영상을 데이터로 표현하는 것은 크게 두 가지 영역으로 나뉘는데, 앞 장까지 했던 화소값이 직접 표현된 공간영역(spatial domain)과 다른 하나가 우주 공간과 같은 변환영역(transform domain)
  • 변환영역은 직교변환에 의해 얻어진 영상 데이터의 다른 표현이다.
    • 여기서는 화소값이 직접 표현되는 것이 아니고 변환계수(coeffcient)로 표현된다.
    • 대표적인 변환은 DCT(Discrete Cosine Transform)와 DFT(Discrete Fourier Transform)이 있는데, 그 중 오래되었고 잘 알려진 것이 DFT인 이산 푸리에 변환이다.
    • 푸리에 변환은 시간(혹은 공간) 영역에서 주파수 영역으로의 변환으로 ‘모든 파형은 단순한 정형파의 합으로 합성되어 질 수 있다’라는 개념에 기초한 해석학적인 방법

공간 주파수의 이해

  • 일반적으로 주파수라는 말은 아래 그림과 같이 1초 동안 진동하는 횟수로 정의한다. 
    • 라디오 방송 채널이나 휴대폰의 통신 대역에서 사용하는 헤르츠(Hz)는 주파수를 표현하는 단위이다.
  • 그러나 이것은 전파라는 신호에 국한된 표현이라 할 수 있다. 아날로그 신호를 디지털화 하는 과정에서 시간단위로 샘플링하는 횟수를 지정할 때에 Hz는 단위와 함께 샘플링 주파수라는 표현을 사용한다.
    • 또한 우리가 공부하는 영상처리에서도 공간 주파수(spatial frequency)라는 개념을 사용한다.

  • 따라서 좀 더 확장된 의미에서 주파수는 이벤트가 주기적으로 발생하는 빈도라고 할 수 있다.
  • 영상에서 화소 밝기의 변화의 정도를 파형의 형태로 그려보면 아래 그림과 같이 표현할 수 있다. 신호의 주파수와 같은 의미가 되는 것이다.
  • 이렇게 확장된 의미를 영상 신호에 적용하면 영상에서의 주파수는 공간상에서 화소 밝기의 변화율이라 할 수 있다.
    • 이런 의미에서 공간 주파수라는 표현을 사용한다.
    • 공간 주파수는 밝기가 얼마나 빨리 변화하는가에 따라 고주파 영역과 저주파 영역으로 분류한다.
  • 아래 그림은 고주파 포함 영역과 저주파 포함 영역을 설명한다.
    • 상단 부분을 보면 화소 밝기가 거의 변화가 없거나 점진적으로 변화하는 것을 볼 수 있는데 이런 부분은 대부분 저주파 성분을 가진 저주파 공간 영역이라 한다.
    • 반면 하단 부분은 화소의 밝기가 급변하는 것을 볼 수 있는데, 이런 부분은 변화가 거의 없는 저주파 성분 위에 변화가 심한 고주파 성분이 포함되어 있는 고주파 공간 영역이라 할 수 있다.
  • 저주파 공간 영역은 보통 영상에서 배경 부분이나 객체의 내부에 많이 있으며, 고주파 공간 영역은 경계부분이나 객체의 모서리 부분에 많이 있다.

  • 변환을 통하여 영상을 주파수 영역별로 분리할 수 있으면, 각 주파수 영역별 처리가 가능할 것이다.
    • 예컨대 경계부분에 많은 고주파 성분을 제거하여 영상을 생성하면 경계가 흐려진 영상을 만들 수 있고, 저주파 성분을 제거하고 고주파 성분만을 취하여 영상을 만든다면 경계나 모서리만 포함하는 영상 즉, 에지 영상이 만들어질 것이다.
  • 일반적인 영상은 공간 영역상에서 저주파 성분과 고주파 성분이 혼합하여 있기 때문에 저주파 영역과 고주파 영역을 분리해서 선별적으로 처리하기란 쉬운 일이 아니다. 따라서 변환영역의 처리가 필요하다.
  • 변환영역, 즉 주파수 영역에서의 영상처리는 아래 그림과 같은 과정을 거친다.
    • 먼저 영상이 입력되면 주파수 영역으로 변환하며,
    • 주파수 변환으로 얻어진 계수의 특정 주파수 영역에 원하는 영상 처리를 적용한다.
    • 마지막으로 처리가 적용된 후에는 다시 주파수 역변환을 통해 공간영역의 영상으로 변환해서 출력 영상을 생성한다.

이산 푸리에 변환

  • 푸리에 변환은 신호나 영상을 주파수 영역으로 변환하는 가장 일반적인 방법으로 다음의 전제를 기본으로 한다.
    • 주기를 가진 신호는 정현파/여현파의 합으로 표현할 수 있다.
  • 여기서 정현파/여현파는 모든 파형 중에 가장 순수한 파형을 말하는 것으로 사인(sin), 코사인(cos) 함수로 된 신호를 말한다. 즉, 사인 또는 코사인 함수의 선형 조합으로 특정 주기의 신호를 구성할 수 있다.
  • 이 전제를 바꾸어 말하면 아래 그림과 같이 주기를 갖는 신호 g(t) 는 여러 개의 사인 및 코사인 함수 g_{f}(t) 로 분리되는 것이다.
    • 여기서 분리된 신호 g_{1}(t), g_{2}(t), g_{3}(t) 는 기저 함수(basis function)가 되며, 기저 함수에 곱해지는 값 0.3, 0.7, -0.5가 주파수의 계수가 된다.
  • 이 계수가 각 주파수 성분의 크기에 해당하며, 신호를 주파수로 변환하는 것은 각 주파수의 기저 함수들에 대한 계수를 찾는 것이다. 또한 주파수 영역에서의 역변환은 각 기저함수와 그 계수들로부터 원본 신호를 재구성하는 것이다.
  • 주파수 변환을 수식으로 표현하면 다음과 같다. 여기서 g_{f}(t) f 주파수에 대한 기저함수이며 G(t) 는 각 기저함수의 계수이다.

g(t) = \int_{-\infty}^{\infty} G(f) \cdot g_{f}(t)df

  • 연속신호에 대한 변환이기 떄문에 -무한대에서 +무한대까지를 적분하여 모든 주파수에 대한 기저함수와 그 계숟르의 선형 조합이 된다.
    • 따라서 연속신호를 주파수 영역으로의 변환은 존재하는 모든 기저함수에 대해 그 계수인 G(t) 를 구하는 것이다.
  • 기저 함수를 어떻게 정하는가에 따라 주파수 영역 변환의 방법이 달라진다.
    • 일반적으로 다음과 같이 사인이나 코사인 함수를 기저함수로 사용한다.이것은 가장 대표적인 방법인 푸리에 변환에서 사용하는 기저함수이다.
    • 또한 푸리에 변환에서는 허수를 이용한 함수를 기저함수로 사용한다.

g_{f}(t) = cos(2 \pi ft) + j \cdot sin(2 \pi ft) = e^{j 2 \pi ft}

  • 이 기저함수로 원본 신호를 나타내면 다음과 같아. 이것은 원본 신호를 만드는 것이기에 푸리에 역변환에 대한 수식이 된다.

g(t) = \int_{-\infty}^{\infty} G(f) \cdot e^{j 2 \pi ft} df

  • 원본 신호로부터 주파수의 계수 G(f) 를 얻는 식은 다음과 같이 유도된다. 이것은 1차원의 연속 신호에 대한 푸리에 변환이다.

G(f) = \int_{-\infty}^{\infty} g(t) \cdot e^{-j 2 \pi ft} dt

  • 이를 디지털 신호에 적용하려면 이산 푸리에 변환(DFT)을 사용해야 한다. 다음은 이산 푸리에 변환과 그 역변환의 수식이다.

G(k) = \sum_{n=0}^{N-1} g[n] \cdot e^{-j 2 \pi k {n \over N}}, (k = 0, ... N-1)

g[n] = {1 \over N} \sum_{n=0}^{N-1} G[k] \cdot e^{j 2 \pi k {n \over N}}, (n = 0, ... N-1)

  • 여기서 g[n] 은 디지털 신호이며, G(k) 는 주파수 k 에 대한 푸리에 변환 계수이다. 
    • 또한 k, n 은 신호의 원소개수(N)만큼 정수로 주어진다.
    • 연속 신호에서 적분기호가 이산신호에서는 합기호로 바뀐다.
  • 2차원 공간상의 영상에 이산 푸리에 변환을 적용하려면 다음 수식과 같이 1차원 이산 푸리에 변환을 가로방향과 세로방향으로 연속해서 두 번 적용해야 한다. 
    • 수식에서 괄호부분이 가로방향에 대한 1차원 푸리에 변환이다.

G(k, l) = \sum_{m=0}^{M-1} (\sum_{n=0}^{N-1} g[n, m] \cdot e^{-j 2 \pi k {n \over N}}) \cdot e^{-j 2 \pi t {m \over M}} \\ = \sum_{m=0}^{M-1} \sum_{n=0}^{N-1} g[n, m] \cdot e^{-j 2 \pi ({kn \over N} + {lm \over M})}

  • 다음 수식은 2차원 이산 푸리에 역변환에 대한 수식이다.

g[n, m] = {1 \over NM} \cdot \sum_{m=0}^{M-1} (\sum_{n=0}^{N-1} G(k, l) \cdot e^{j 2 \pi k {n \over N}}) \cdot e^{j 2 \pi l {m \over M}} \\ = {1 \over NM} \sum_{m=0}^{M-1} \sum_{n=0}^{N-1} G(k, l) \cdot e^{j 2 \pi ({kn \over N} + {lm \over M})}

  • 푸리에 변환의 기저함수에 허수를 이용했기 때문에 실수부와 함께 허수부에 대한 고려도 해야 한다.
    • 다음은 기저함수를 사인과 코사인 함수로 변경하여 1차원 푸리에 변환을 수식으로 나타내면 다음과 같다.

G(k) = \sum_{n=0}^{N-1} g[n] \cdot (cos(-2 \pi k {n \over N}) + j \cdot sin(-2 \pi k {n \over N}))

  • 그리고 실수부와 허수부를 구분하여 표현하면 다음과 같다. 푸리에 변환과 그 역변환이 사인과 코사인 함수에서 각도의 부호만 반대이다.

G(k)_{Re} = \sum_{n=0}^{N-1} g[n]_{Re} \cdot cos(-2 \pi k {n \over N}) - g[n]_{Im} \cdot sin(-2 \pi k {n \over N})

G(k)_{Im} = \sum_{n=0}^{N-1} g[n]_{Im} \cdot cos(-2 \pi k {n \over N}) + g[n]_{Re} \cdot sin(-2 \pi k {n \over N})

g[n]_{Re} = {1 \over N} \sum_{n=0}^{N-1} G(k)_{Re} \cdot cos(2 \pi k {n \over N}) - G(k)_{Im} \cdot sin(2 \pi k {n \over N})

g[n]_{Im} = {1 \over N} \sum_{n=0}^{N-1} G(k)_{Im} \cdot cos(2 \pi k {n \over N}) + G(k)_{Re} \cdot sin(2 \pi k {n \over N})

  • 푸리에 변환을 수행하면 복소수의 행렬이 결과로 생성된다. 이것을 영상으로 확인하기 위해서는 복수부의 실수부와 허수부를 벡터로 간주하여 다음의 수식과 같이 벡터의 크기를 구하면 된다. 이것을 주파수 스펙트럼이라 한다.
    • 또한 실수부와 허수부의 각도를 이용해서 주파수 위상을 계산할 수도 있다.

|G(k, l)| = \sqrt{Re(k, l)^{2} + Im(k, l)^{2}}

\theta (k, l) = tan^{-1} [ {Im(k, l) \over Re (k, l)} ]

  • 여기서 주파수 스펙트럼 영상은 저주파 영역의 계수값이 고주파 영역에 비해 상대적으로 너무 크다.
    • 이로 인해 계수값을 일반적인 방법으로 정규화해서 영상으로 표현하면 최저 주파 영역만 흰색으로 나타나고 나머지 영역은 거의 검은색으로 나타나서 고주파 영역의 계수를 영상으로 확인하기가 곤란하다.
    • 이런 문제를 해결하기 위해 계수값에 로그 함수를 먼저 적용하고 정규화한다.
void log_mag(Mat complex_mat, Mat& dst)
{
Mat planes[2];
split(complex_mat, planes); // 2채널 복소행렬 분리
magnitude(planes[0], planes[1], dst); // 크기 계산
log(dst + 1, dst);
normalize(dst, dst, 0, 255, CV_MINMAX); // 정규화 수행
dst.convertTo(dst, CV_8U);
}
  • 위의 log_mag() 함수는 DFT를 수행한 행렬에서 주파수 스펙트럼을 계산하고, 로그 함수를 적용한다.
    • 여기서 첫 번째 인수인 DFT 결과 행렬(complex_mat)은 실수부와 허수부를 갖는 2채널 행렬이다.
    • 따라서 cv::split() 함수로 2채널 행렬을 1채널 행렬 2개로 분리해서 크기를 구한다.
  • 또한 DFT 수행 후의 주파수 스펙트럼 영상은 저주파 영역이 영상의 모서리 부분에 위치하고 고주파 부분이 중심부에 있다.
    • 즉, 아래 그림의 왼쪽과 같이 푸리에 변환 영상에서 사각형의 각 모서리를 중심으로 원형의 밴드를 형성하여 주파수 영역이 분포한다. 
    • 이 때문에 해당 주파수 영역에서 어떤 처리를 하려면 상당한 불편함이 있다. (여기서 원형의 밴드는 이해를 돕기 위해 표현한 것으로 실제 DFT 스펙트럼 영상에서 나타나는 것은 아니다.)
  • 이 문제는 1사분면과 3사분면의 영상을 맞바꾸고, 2사분면과 4사분면의 영상을 맞바꿈으로서 해결할 수 있다.
    • 결과적으로 아래 그림의 오른쪽과 같이 영상의 중심이 최저주파 영역, 그리고 바깥쪽이 고주파 영역이 되며, 그림과 같이 원형의 밴드로 주파수 영역을 쉽게 구분할 수 있다.
    • 이런 과정을 셔플링(shuffling) 혹은 시프트(shift) 연산이라고 한다.
  • 다음의 shuffling() 함수는 입력 행렬에 셔플링을 수행해서 반환하는 함수이다.
    • q1~q4의 Rect 객체를 통해서 각 사분면의 영역을 지정한다.
    • 그리고 원본 행렬(mag_img)과 반환 행렬(dst) 모두에서 각 사분면을 관심영역으로 지정하여 참조한다.
    • Mat::copyTo() 함수를 통해서 참조된 원본 행렬을 참조된 반환 행렬로 복사하여 사분면의 맞교환을 수행한다.
void shuffling(Mat mag_img, Mat& dst)
{
int cx = mag_img.cols / 2;
int cy = mag_img.rows / 2;

Rect q1(cx, 0, cx, cy); // 1사분면 사각형
Rect q2(0, 0, cx, cy); // 2사분면 사각형
Rect q3(0, cy, cx, cy); // 3사분면 사각형
Rect q4(cx, cy, cx, cy); // 4사분면 사각형

dst = Mat(mag_img.size(), mag_img.type());
mag_img(q1).copyTo(dst(q3));
mag_img(q3).copyTo(dst(q1));
mag_img(q2).copyTo(dst(q4));
mag_img(q4).copyTo(dst(q2));
}
  • N x M 크기의 영상에서 2차원 푸리에 변환은 N2 x M2 만큼의 시간 복잡도를 요구한다.
    • 따라서 영상의 크기가 커지면 수행속도는 기하급수적으로 증가한다.
    • 이러한 이유에서 푸리에 변환을 빠르게 수행하는 알고리즘의 필요성이 대두되었다.

고속 푸리에 변환(FFT: Fast Fourier Transform)

  • 이산 푸리에 변환은 원본 신호의 한 원소에 곱해지는 기저 함수의 원소들을 원소 길이만큼 반복적으로 곱해야하기 때문에 신호가 커질수록 계산 속도는 기하급수적으로 증가한다.
    • 고속 푸리에 변환은 이 과정을 삼각함수의 주기성을 이용해 작은 단위로 분리해서 반복적으로 수행하고 합치도록 하여 효율성을 높이는 방법이다.
  • 다음은 삼각함수의 주기성을 이용하는 방법을 간단히 설명한다. 먼저 푸리에 변환 수식에서 짝수 번째 부분(2n)과 홀수 번쨰 부분(2n+1)을 분리하여 수식을 다음과 같이 정리한다.

G(k) = \sum_{n=0}^{L-1} g[2n] \cdot e^{-j 2 \pi k {2n \over 2L}} + \sum_{n=0}^{L-1} g[2n+1] \cdot e^{-j 2 \pi k {2n + 1 \over 2L}}

G(k) = \{ \sum_{n=0}^{L-1} g[2n] \cdot e^{-j 2 \pi k {n \over L}} \} + \{ \sum_{n=0}^{L-1} g[2n+1] \cdot e^{-j 2 \pi k {n \over L}} \} \cdot e^{-j 2 \pi k {1 \over 2L}}

  • 여기서 짝수 신호와 홀수 신호를 다음과 같이 지정해 보자

G_{even}(k) = \sum_{n=0}^{L-1} g[2n] \cdot e^{-j 2 \pi k {n \over L}}

G_{odd}(k) = \sum_{n=0}^{L-1} g[2n+1] \cdot e^{-j 2 \pi k {n \over L}} 

  • 그러면 푸리에 변환 수식은 다음과 같이 나타낼 수 있다.

G(k) = G_{even}(k) + G_{odd}(k) \cdot e^{-j 2 \pi k {1 \over 2L}}

  • 위 수식의 공통 지수부분에서 한 주기(L)를 더한 수식을 정리해 보자. 이것은 삼각 함수의 주기성으로 인해 다음과 같이 뒷부분 지수를 제거할 수 있다.

e^{-j 2 \pi (k+L) {n \over L}} \\ = e^{-j 2 \pi k {n \over L}} \cdot e^{-j 2 \pi  {Ln \over L}} \\ = e^{-j 2 \pi k {n \over L}} \cdot e^{-j 2 \pi n} \\ = e^{-j 2 \pi k {n \over L}}

  • 이 주기성을 이용하여 G(k+L)을 계산하면 다음과 같다.

G(k+L) = G_{even}(k+L) + G_{odd}(k+L) \cdot e^{-j 2 \pi (k+L) {1 \over 2L}} \\ = G_{even}(k) + G_{odd}(k) \cdot e^{-j 2 \pi (k+L) {1 \over 2L}}

  • 여기서 e^{-j 2 \pi (k+L) {1 \over 2L}} 의 지수를 분리하여 정리하면 다음과 같다.

e^{-j 2 \pi (k+L) {1 \over 2L}} \\ = e^{-j 2 \pi k {1 \over 2L}} \cdot e^{-j 2 \pi  {L \over 2L}} \\ = e^{-j 2 \pi k {1 \over 2L}} \cdot e^{-j \pi} \\ = -e^{-j 2 \pi k {1 \over 2L}}

  • 따라서 최종적으로 G(k+L)은 다음과 같다.

G(k+L) = G_{even}(k) - G_{odd}(k) \cdot e^{-j 2 \pi k {1 \over 2L}}

  • 이것은 G(k+L)의 값이 G(k)의 값들을 이용해서 구할 수 있다는 것이다.
  • 아래 그림은 삼각함수의 주기성이 어떻게 적용되어 반복계산을 줄이는지를 보인다.
    • 8개의 원소를 갖는 원본 신호 g[n]을 짝수 신호와 홀수 신호로 구분하여 푸리엔 변환을 한다.
    • 짝수 DFT 결과와 홀수 DFT 결과로 최종 변환 신호 G[n]을 구성하는 방법을 보인다.
  • 짝수 원소와 홀수 원소를 한 번 분리하는 것으로는 수행속도를 줄이지 못한다.
    • 각 그룹 내 원소들을 연속적으로 분리할 수 있을 것이다. 즉 아래 그림과 같이 짝수 원소 그룹내에서 짝수 원소와 홀수 원소를 다시 분리하고, 홀수 원소 그룹 내에서 짝수 원소와 홀수 원소를 분리한다.
    • 이렇게 연속적으로 분리하면 최종적으로 입력 신호를 2개 원소씩 묶을 수 있다.
  • 여기서 입력 신호에 대해 짝수부와 홀수부로 계속적으로 분리하여 최종적으로 두 원소만 갖게끔 신호를 재배열해야 한다. 이것을 스크램블(scramble)이라고 한다.
    • 스크램블은 보통 비트의 순서를 바꾸는 방법을 설명하지만 속도면에서 비효율적이며, 다음의 scramble() 함수로 구현한다.
Mat scramble(Mat signal)
{
Mat dst = signal.clone();

for (int i = 0, j = 0; i < dst.cols - 1; i++)
{
if (i > j)
{
swap(dst.at<Vec2f>(i), dst.at<Vec2f>(j));
}

int m = dst.cols >> 1;

while ((j >= m) && (m >= 2))
{
j -= m;
m >>= 1;
}

j += m;
}

return dst;
}
  • 다음은 버터플라이(butterfly) 과정이다. 이것은 스크램블 결과 원소에서 이웃한 두 원소에 대해 이산 푸리에 변환을 수행하는 것이다.
    • 여기서 Wk는 푸리에 변환의 기저함수인 삼각함수의 수식이다.
    • 버터플라이는 아래 그림과 같이 흐름도의 모양이 나비와 비슷해서 붙여진 이름이다.
  • 버터플라이 과정은 원본 신호 길이를 두 개 원소 신호로 분리하며, 분리 횟수만큼 연속적으로 반복한다.
    • 여기서 원본 신호를 연속적으로 짝수부와 홀수부로 분리하기 때문에 원본 신호의 원소 개수는 2의 자승이 되어야 한다.
    • 영상의 크기가 반드시 2의 자승이 되는 것이 아니기 때문에 이 문제를 해결하는 방법으로 원본 영상의 가로와 세로 크기를 2의 자승이 되게 넓히고, 빈 공간을 검은색(0)으로 채우는 방법을 사용한다. 이것을 영삽입(zero-padding)이라 한다.
  • 영삽입이 가능한 것은 원본 영상의 평행이동이 푸리에 변환 결과에 영향을 미치지 않기 때문이다.
    • 즉, 아래 그림의 오른쪽 영상에서 원본 영상이 어느 위치로 평행이동되든지 상관없이 푸리에 변환 결과에서 스펙트럼은 동일하다.
  • 다음의 zeropadding() 함수는 입력 영상에 영삽입을 수행해서 반환하는 함수이다. 그 과정은 2에 대한 로그함수인 log2() 함수로 영상의 가로와 세로가 2의 몇 승인지를 계산한다.
    • 계산된 승수는 소수점을 포함하며, ceil() 함수로 소수점 부분을 올림처리한다.
    • 그리고 다시 계산된 승수에 << 연산으로 2의 자승을 만든다.
    • 이렇게 하면 원본 영상보다 소수점 올림난큼 큰 2의 자승 크기(m, n)을 계산할 수 있다.
Mat zeropadding(Mat img)
{
int m = 1 << (int)ceil(log2(img.rows));
int n = 1 << (int)ceil(log2(img.cols));
Mat dst(m, n, img.type(), Scalar(0));

Rect rect(Point(0, 0), img.size());
img.copyTo(dst(rect));
dst.convertTo(dst, CV_32F);

return dst;
}

FFT를 이용한 주파수 영역 필터링

주파수 영역 필터링의 과정

  • 영상을 주파수 영역으로 변환하면 화소의 밝기가 서서히 변화하는 저주파 영역과 급격하게 변화하는 고주파 영역을 공간 영역에 비하여 쉽게 분리할 수 있다.
    • 이렇게 분리된 주파수 영역에 대해 각 주파수 영역을 강화하거나 약화하거나 혹은 제거하는 등의 처리를 통해 다양한 영상 처리를 할 수 있다.
  • 앞서 배운 내용으로 영상에 2차원 푸리에 변환을 수행할 수 있다. 이를 통해 아래 그림과 같이 영상을 주파수 영역으로 쉽게 변환할 수 있다.
    • 주파수 영역에서 필터링 과정은 푸리에 변환 계수에 필터 행렬을 원소간(element-wise)에 곱하여 수행된다.
    • 여기서 푸리에 변환 계수는 복소수이기 때문에 필터의 곱셈도 실수부와 허수부의 두 채널에 수행해야 한다.
    • 마지막으로 필터링된 푸리에 변환 계수를 푸리에 역변환(IFFT)함으로써 다시 공간영역의 영상으로 만들 수 있다.
  • 이러한 일련의 과정을 주파수 성분 조작이라 한다. 또한 필터를 어떻게 구성하느냐에 따라 저주파 통과 필터링, 고주파 통과 필터링, 대역 통과 필터링 등을 쉽게 구현할 수 있다.

저주파 및 고주파 통과 필터링

  • 저주파 통과 필터링은 DFT 변환 영역에서 저주파 영역의 계수들은 통과시키고, 그 외의 영역 즉, 고주파 영역의 계수는 차단하는 것을 말한다.
    • 아래 그림의 왼쪽처럼 푸리에 변환을 하고, 셔플링을 수행해서 주파수 스펙트럼 영상을 보면 중심 부분이 저주파 영역이며, 외곽으로 갈수록 고주파 영역이다. 여기서 원형의 점선은 주파수 밴드를 시각적으로 표현한 것이다.
    • 필터링은 주파수 계수에 필터 행렬의 원소가 곱해져서 수행된다. 따라서 저주파 통과 필터의 모양은 아래 그림의 중간 그림처럼 중심에서 지정된 반지름만큼 원형으로 1의 값을 갖게 하고, 외곽 부분을 0으로 지정하면 된다. 그림에서 흰색은 1의 값이며, 검은색은 0의 값이다.
    • 고주파 통과 필터링은 저주파 통과 필터와는 반대로 고주파 영역의 계수들을 통과시키고 저주파 영역의 계수들은 차단하는 것이다. 아래 그림의 오른쪽과 같이 중심에서 지정된 반지름 크기의 원형으로 0의 값을 갖게 하고 가장자리 부분을 1로 지정하면 된다.
    • 주파수 계수와 필터의 원소가 곱해지기 때문에 0의 값을 갖는 부분은 주파수 계수가 제거되어 차단되고, 1의 값을 갖는 부분은 그대로 유지되어 통과된다.

버터워스, 가우시안 필터링

  • 앞선 대역 통과 필터는 특정한 대역에서 급격하게 값을 제거하기 때문에 결과 영상의 화질이 좋지 못하다.
    • 특히 저주파 통과 필터링의 경우 영상에서 객체의 경계부분이 완만해지기는 하지만, 경계부분 주위로 잔물결 같은 무늬 (ringing pattern)가 나타나서 화질이 더욱 떨어진다.
    • 이 문제를 해결하는 것은 필터 원소의 값을 차단 주파수에서 급격하게 0으로 만들지 않고 완만한 경사를 이루도록 구성하면 된다.
    • 대표적으로 버터워즈 필터링(Butterworth filter)과 가우시안 필터(Gaussian filter)가 있다.
  • 가우시안 필터는 필터 우너소의 구성을 가우시안 함수의 수식 분포를 갖게 함으로써 차단 주파수 부분을 점진적으로 구성한 것이다.
    • 가우시안 함수의 수식에서 표준편차(\sigma )를 주파수를 차단할 반지름의 위치(R )로 간주하자.
    • 그러면 아래 그림의 오른쪽과 같이 포물선의 곡선을 갖는 주파수 공간 필터가 구성된다.
    • 또한 가우시안 함수의 수식은 {1 \over 2 \pi \sigma^{2}} 를 곱해야 하지만, 원소 값의 스케일을 최댓값이 1이 되도록 하기 위해 생략한다.

f(x, y) = exp(-{dx^{2} + dy^{2} \over 2 R^{2}}) \\ dx = x - center x \\ dy = y - center y

  • 버터워즈 필터는 다음의 수식으로 필터 원소의 구성이 가능하다. 여기서 차단 주파수 반지름 위치(R )와 지수의 승수인 n 값을 어떻게 지정하느냐에 따라 차잔 필터의 반지름과 포물선의 곡률이 달라진다.

f(x, y) = -{1 \over 1 + ({\sqrt{dx^{2} + dy^{2}} \over R})^{2n}} \\ dx = x - center x \\ dy = y - center y

이산 코사인 변환

  • 1974년 미국 텍사스 대학에서 라오 교수팀이 이산 코사인 변환(DCT: Discrete Cosine Transform)이라는 새로운 직교변환에 관한 논문을 발표하면서 멀티미디어 혁명이 시작되었다.
    • 라오 팀은 영상 신호의 에너지 집중 특서잉 뛰어나서 영상 압축에 효과적인 주파수 변환 방법을 찾는 것이 목표였고 그 결과가 바로 DCT인 것이다.
  • 이산 푸리에 변환(DFT)은 실수부에 코사인 함수가 곱해지며, 허수부에 사인 함수가 곱해져서 이루어진다. 반면 이산 코사인 변환은 이산 푸리에 변환에서 실수부만 취하고, 허수부분을 제외함으로써 코사인 함수만으로 구성된 직교 방법이다. 이를 이산 여현변환이라고도 한다.
  • 다음은 1차원 이산 코사인 변환의 수식이다.

F(k) = C(k) \cdot \sum_{n=0}^{N-1} g[n] \cdot cos ({(2n + 1) k \pi \over 2N})

g[n] = \sum_{k=0}^{N-1} C(k) \cdot F(k) \cdot cos ({(2n + 1) k \pi \over 2N})

$latex k = 0, 1, … N-1 \\ C(k) \Rightarrow k = 0 \Rightarrow \sqrt{1 / N}, k \neq 0 \Rightarrow \sqrt{2 / N} &s=2$

  • F(k) 는 주파수 영역 신호이며, g(k) 은 공간 영역의 신호이다.

F(k, l) = C(k) \cdot C(l) \cdot \sum_{n=0}^{N-1} \sum_{m=0}^{m-1} g[n, m] \cdot cos ({(2n + 1) k \pi \over 2N}) \cdot cos ({(2m + 1) l \pi \over 2M})

g(n, m) = \sum_{k=0}^{N-1} \sum_{l=0}^{m-1} C(k) \cdot C(l) \cdot F(k, l) \cdot cos ({(2n + 1) k \pi \over 2N}) \cdot cos ({(2m + 1) l \pi \over 2M})

k = 0, 1, ... N-1 \\ C(k) \Rightarrow k = 0 \Rightarrow \sqrt{1 / N}, k \neq 0 \Rightarrow \sqrt{2 / N}

l = 0, 1, ... , M-1 \\ C(l) \Rightarrow l = 0 \Rightarrow \sqrt{1 / N}, l \neq 0 \Rightarrow \sqrt{2 / N}

  • DCT는 일반적으로 전체 영상을 한 번에 변환시키는 것이 아니라 영상을 작은 블록으로 나누어서 블록 단위로 수행한다.
    • 이 블록의 크기를 키울수록 압축의 효율이 높아지지만, 변환의 구현이 어려워지고 속도도 느려진다.
    • 일반적으로 8 x 8 크기가 성능과 구현 용이성간의 상호보환(trade-off)되어 표준으로 사용된다.
  • 8 x 8 크기의 한 블록을 DCT 변환하면 아래 그림과 같이 64개의 주파수 계수가 구성된다. (0, 0) 위치를 DC 계수라 하며, 나머지 계수들을 AC 계수라 한다.
    • DC 계수는 공간 영역의 화소값의 평균에 해당하는 값으로서 에너지가 집중되어 있고, 영상의 주요 성분을 포함하고 있다.
    • 각 주파수 계수는 영상의 밝기 변화 특성을 나타낸다. 왼쪽 상단으로 갈수록 저주파 영역이며, 오른쪽 하단으로 갈수록 고주파 영역이다. 저주파 영역으로 갈수록 밝기 변화가 적으며, 고주파 영역으로 갈수록 밝기 변화의 정도가 증가한다.
  • 블록의 크기를 N x M이라 할 때, 코사인 함수 부분이 블록의 한 화소에서 N x M 만큼 계산하고, 블록의 모든 화소에서 코사인 함수를 계산하게 되면 O(N2 x M2)의 계산 복잡도가 된다.
    • 그리고 원본 영상을 N x M 크기의 블록으로 나누어 모든 블록에 DCT를 수행한다. 따라서 DCT 변환의 전체 계산복잡도는 다음의 수식과 같다.

DCT 변환의 계산 복잡도 = $latex O(N^{2} \times M^{2} \times Total_{block}) \\ Total_{block} = H_{block} \times W_{block} \\ H_{block} = {ImageHeight \over N} \\ W_{block} = {ImageWidth \over M} &s=2$

  • 512 x 512 영상에 8 x 8 블록으로 DCT 를 수행하면, 한 블록의 DCT 복잡도가 4,096번이며, 전체 블록의 개수가 4,096개이므로 16,777,216번의 코사인 함수를 계산해야 한다. 게다가 코사인 함수의 계산이 상대적으로 느리기 때문에 상당한 시간이 소요된다.
    • 다행히 DCT 정변환과 역변환에서 코사인 함수 부분이 동일하고 한 블록에 수행되는 모든 코사인 함수 값은 블록마다 반복적으로 사용되므로 블록에 사용되는 코사인 함수값들을 미리 계산에서 행렬에 저장해 두면 효율적이다.
    • (관련 코드 생략)
[ssba]

The author

지성을 추구하는 사람/ suyeongpark@abyne.com

댓글 남기기

This site uses Akismet to reduce spam. Learn how your comment data is processed.