주파수 영역에서의 영상처리

1. 개요
영상의 공간적인 분석이 아니라 영상을 이루고 있는 주파수 영역으로의 변환을 통해서도 영상내의 또 다른 유용한 정보를 추출할 수 있다. 아래 그림은 다양한 영상에 대한 이산 푸리에 변환을 통한 주파수 영역으로의 영상의 변환 예이다. 원 영상은 위의 열에 해당하며 영상 하나하나가 주파수 영역으로의 변환 결과는 바로 밑 열에 나타나 있다.
위 그림의 결과를 보면 원 영상에서의 복잡한 밝기값의 변화 정도가 주파수 영역으로의 변환 후에는 현저히 줄어듦을 알 수 있다. 위의 세 영상 모두 저주파 성분(변환 후의 중앙 부분)이 지배적임을 알 수 있다. 이러한 성질을 이용한 것이 최근에 영상 압축 표준 등에서 널리 쓰이고 있는 DCT(Discrete Consine Transformation)이다.

2. 이론적인 배경
영상 신호의 경우 삼차원 공간상에 발생하는 사건을 2차원 평면상의 특정 시간에 매핑한 것으로 불 수 있으며, 이해의 편의를 돕기 위해 이를 1차원 신호에 대해 생각하기로 한다.

위의 식은 연속 신호인 경우에 해당하며, 이를 이산 신호에 적용하기 위해서는 먼저 연속 신호를 시간축에 대해 샘플링할 필요가 있다. 이와 함께 각 시간상에서의 신호의 크기도 컴퓨터상에서 사용할 수 있도록 일정 값을 가져야 한다(양자화). 연속 신호를 이산 신호로 변경한 후 이산 푸리에 변환(DFT : Discrete Fourier Transformation)을 이용하여 신호의 주파수 영역에서의 특성을 알아낼 수 있다.

N개의 신간축 상에서 샘플링한 이산 신호값이 주어진 경우 이산 푸리에 변환을 통해 N개의 주파수 상에서의 변환값을 알 수 있다. 이산 푸리에 변환은 다음의 수식과 같다.
주파수 영역에서 시간 영역으로의 변환은 다음 식과 같다.
위의 두 식을 이용하면 시간축 상의 이산 신호를 주파수 영역으로 변경하고 또한 그 역과정도 계산할 수 있다. 식 9-2를 이용하여 계산할 경우, S(0),…,S(N-1) 각각의 계산시 N개의 곱셈과 (N-1)번의 덧셈이 필요함을 알 수 있다. 이로부터 N개의 점이 주어진 경우 이산 푸리에 변환의 경우 N제곰의 곱셈이 필요함을 알 수 있다.3. 고속 푸리에 변환
이 방법은 식 9-2의 실제 계산에서 나타나는 항의 주기성과 대칭성을 이용하고 있다. 주어진 신호의 개수 N이 다음을 따른다고 가정한다.
N은 짝수개이므로 이를 이등분할 수 있으며 이는 다음과 같이 나타낼 수 있다.
식 9-4를 식 9-2에 대입하면 다음과 같다.
식 9-5에서 식 9-5에서 N점의 계산시 이를 각각 짝수 번째와 홀수 번째의 그룹으로 나누어 계산할 수 있음을 알 수 있다.N점의 계산시 이를 각각 짝수 번째와 홀수 번째의 그룹으로 나누어 계산할 수 있음을 알 수 있다.
식 9-6을 식 9-5에 대입하면 다음과 같다.
식 9-7을 구성하고 있는 두 항목 각각에 대해 다음과 같이 정의할 수 있다.
식 9-8에서 항은 주어진 N개의 점들 중 짝수 번째에 해당하는 점둘 N/2개를 이용하여 계산한 값이며, 항은 주어진 N개의 점들 중 짝수 번째에 해당하는 점들 N/2개를 이용하여 계산한 값이다.

식 9-8을 식 9-7에 대입하면 다음과 같다
식 9-9를 이용하여 N개의 점이 주어진 경우 이를 각각 짝수 번째 점들과 홀수 번째 점들의 두 그룹으로 나누어 계산과정을 수행 후 S(0),…,S(L-1) = S(N/2-1)까지의 값을 구할 수 있다. 나머지 S(L)=S(N/2),S(L+1),…,S(N-1)의 값은 다음의 식으로부터 구할 수 있다.
식 9-10을 통해 나머지 절반의 점들의 계산시 기존의 계산 결과를 그대로 이용할 수 있음을 알 수 있으며, 이를 통해 고속 연산이 가능하게 된다. 위의 경우에는 N개의 신호점들을 N/2개의 짝수 번째 신호점들과 N/2개 홀수 번째 그룹으로 나눈 것이며, 동일하게 N/2개의 짝수 번째 신호점들을 또 다시 N/4개의 짝수 그룹과 N/4개의 홀수 그룹으로 나눌 수 있다. 이와같은 과정을 최적적으로 2개의 신호점으로 이루어진 그룹까지 수행 후 식 9-9와 식 9-10을 이용하여 계산하게 된다.

8개의 신호점이 주어진 경우에 대해 예를 들어보기로 한다. 아래 그림에서 트리의 왼쪽 자식은 짝수번째 신호점들의 그룹이며 올느쪽 자식은 홀수 번째 신호점들의 그룹이다. 최종적으로 2개의 신호점으로 이루어진 그룹까지 트리가 분기함을 알 수 있다. 그림에서 보는 바와 같이 식 9-9와 식 9-10을 이용한 고속 푸리에 계산시 입력 신호의 재배치 과정이 필요함을 알 수 있다. 트리의 최하단에서 각각 2점을 이용한 계산을 수행 후 이를 트리의 상단으로 전파시키면서 각각 4점에 의한 연산, 8점에 의한 연산 등의 결과를 얻을 수 있다. 이를 통해 기존의 계산 결과를 이용함으로써 고속 연산이 가능하게 된다.

아래 그림의 비교에서 알 수 있듯이 8점에 대한 이산 푸리에 변환의 계산시 식 9-2의 단순 적용에 의한 계산 결과인 그림 (a)는 8 * 8 = 64번의곱셈이 필요하며, 식 9-9와 식 9-10에 의한 고속 푸리에 계산의 경우 (b)에서와 같이 8 * 3 = 24번의 곱셈이 필요하다.
지금까지는 1차원의 이산 신호에 대해서 다루었으며, 우리가 다루는 영상 이미지의 경우 1차원 신호가 아니라 2차원 신호의 형태이므로 주파수 영역으로의 변환시 2차원 함수에 대한 이산 푸리에 변환(DFT)을 이용한다.
식 9-11은 2차원 이산 푸리에 변환(DFT)을 나타내고 있으며, 위의 식 중 2번째 열에서 나타난 바와 같이 2차원 이산 푸리에 변환의 계산시에는 분리 성질을 이용할 수 있음을 알 수 있다. 먼저 열이나 행 중 하나의 방향으로 고속 푸리에 변환(FFT)을 수행하면 2차원 이산 푸리에 변환(DCF) 결과를 얻을 수 있다. 2차원 이산 푸리에 변환의 계산시 이와 같은 성질을 이용하여 1차원 고속 푸리에 변환을 두 번 수행하면 된다. 1차원 신호에 대한 고속 푸리에 변환에 관한 코드는 아래 리스트와 같다. 리스트의 전반부는 그림 9-2와 같은 입력 신호의 재배치를 수행하는 과정이며 후반부는 재배치된 신호를 이용하여 고속 푸리에 변환을 구현하고 있다.

/*NumData 2의 지수승이어야 한다 */ 
/*Foward: 1인 경우 Discrete Fourier Transform (Forward) */ 
/*        0인 경우 Inverse Discrete Fourier Transform (Backward) */ 
void CWinTestDoc::m_FFT1D(int NumData, float *pOneDRealImage, float *pOneDImImage, int Forward)
{
  int log2NumData;
  int HalfNumData;
  int i,j,k,mm;
  int step;
  int ButterFly,ButterFlyHalf;
  float RealValue,ImValue;
  double AngleRadian;
  float CosineValue,SineValue;
  float ArcRe,ArcIm,ReBuf,ImBuf,ArcBuf;

  /*scramble 과정 수행 */ 
  /*입력 데이터의 순서를 바꿈 */ 
  log2NumData=0;
  while(NumData != (1<<log2NumData))
    log2NumData++;
  HalfNumData=NumData>>1;
  j=1;
  for(i=1;i<NumData;i++)
  {
    if(i<j)
    {
      RealValue=pOneDRealImage[j-1];
      ImValue=pOneDImImage[j-1];
      pOneDRealImage[j-1]=pOneDRealImage[i-1];
      pOneDImImage[j-1]=pOneDImImage[i-1];
      pOneDRealImage[i-1]=RealValue;
      pOneDImImage[i-1]=ImValue;
    }
    k=HalfNumData;
    while(k<j)
    {
      j -= k;
      k=k>>1;
    }
    j += k;
  }
  /*butterfly 과정 수행 */ 
  for(step=1;step<=log2NumData;step++)
  {
    ButterFly=1<<step;
    ButterFlyHalf=ButterFly>>1;
    ArcRe=1;
    ArcIm=0;
    AngleRadian=(double)(3.141592/ButterFlyHalf);
    CosineValue=(float)cos(AngleRadian);
    SineValue=(float)sin(AngleRadian);
    if(Forward) /*Foward*/ 
      SineValue=-SineValue;
    for(j=1;j<=ButterFlyHalf;j++)
    {
      i=j;
      while(i<=NumData)
      {
        mm=i+ButterFlyHalf;
        ReBuf=pOneDRealImage[mm-1]*ArcRe-pOneDImImage[mm-1]*ArcIm;
        ImBuf=pOneDRealImage[mm-1]*ArcIm+pOneDImImage[mm-1]*ArcRe;
        pOneDRealImage[mm-1]=pOneDRealImage[i-1]-ReBuf;
        pOneDImImage[mm-1]=pOneDImImage[i-1]-ImBuf;
        pOneDRealImage[i-1]=pOneDRealImage[i-1]+ReBuf;
        pOneDImImage[i-1]=pOneDImImage[i-1]+ImBuf;
        i += ButterFly;
      }
      ArcBuf=ArcRe;
      ArcRe=ArcRe*CosineValue-ArcIm*SineValue;
      ArcIm=ArcBuf*SineValue+ArcIm*CosineValue;
    }
  }
  if(Forward) /*Forward*/ 
  {
    for(j=0;j<NumData;j++)
    {
      pOneDRealImage[j] /= NumData;
      pOneDImImage[j] /= NumData;
    }
  }
}