해당 포스트는 "열혈강의 영상처리 프로그래밍" 책의 내용을 요약한 것이다.
해당 포스트 읽기 전에 이전 포스트를 읽기를 추천한다.
※ 고속 푸리에 변환
고속 푸리에 변환은 이산 푸리에 변환의 기본적인 구현에서 반복해서 수행하는 계산을 한 번씩만 수행하도록 하여 효율성을 높이는 방법이다. 이산 푸리에 변환을 계산할 때 사인과 코사인 함수를 곱하는데 사인과 코사인 함수는 주기성을 가져 다른 정의역을 대입해도 같은 결과가 나올 수 있다. 그래서 이산 푸리에 변환을 할 때 기저함수들의 같은 값들이 계속 곱해질 수 있다. 이를 개선한 것이 고속 푸리에 변환이다.
고속 푸리에 변환은 신호의 길이 N이 N=2^d와 같이 2의 거듭제곱 형태라고 가정한다. N은 짝수이므로 N=2L과 같이 나타낼 수 있고 이를 이산 푸리에 변환 식에 적용하면 다음과 같다.
위 식에서 n이 짝수인 경우와 홀수인 경우로 나누면 다음과 같다.
위 식에서 마지막 등호 대괄호 두 개 있는 부분에서 왼쪽 대괄호는 입력 신호 g[2n]에 대한 독립적인 푸리에 변환에 해당하고 오른쪽 대괄호는 g[2n+1]에 대한 독립적인 푸리에 변환에 해당한다. 그러면 다음과 같이 짝수 신호와 홀수 신호로 표현 할 수 있다.
위 두 표현을 사용해서 원래의 G(k)를 다음과 같이 나타낼 수 있다.
고속 푸리에 변환은 이산 푸리에 변환에서 삼각함수의 주기성을 이용해서 더 효율적인 계산을 이루도록 한다고 했다. 기저함수 
위 식에서 1인 자연 지수 값은 푸리에 변환에서 기저 함수 사인,코사인 함수를 자연 지수 함수로 바꾼 식을 이용하여 자연 지수를 사인,코사인 함수로 바꾸면 1이 나오게 된다. 위 식을 이용하면 다음 식이 성립하게 된다.
이를 이용해 G(K+L) 값을 구해보면 다음과 같다.
위 식의 자연 지수에 k = k+L 을 대입하면 다음과 같다.
위에서 자연 지수를 삼각함수의 주기성에 의해 변환한 것과 같이 위와 같은 등식이 나올 수 있다. 위 등식을 G(k+L)에 대입하면 다음과 같은 식이 나온다.
위 식이 의미하는 것은 주파수 함수 G(k)의 k번째 원소를 구할 때 썼던 값들을 이용하여 G(k+L) 값을 바로 구할 수 있다는 것이다. 이산 푸리에 변환에서는 전체 주파수 함수값을 구하는 계산량은 입력 신호의 원소 수 N의 제곱에 비례하는 데 이렇게 짝수와 홀수 원소로 분류하여 변환을 수행하면 원소 수가 N/2인 신호를 두 번 변환하는 것이 되어 계산량이 2*(N/2)*(N/2)에 비례하게 된다. 아래는 짝수 신호와 홀수 신호로 분류한 이산 푸리에 변환을 길이가 8인 입력 신호 f(n)에 대해 수행하는 모습이다.
그림2.
보다시피 짝수와 홀수로 분류해서 이산 푸리에 변환을 수행하고 있고 F(4)와 F(0)의 값을 계산하는 데 똑같이 F<even>(0)과 F(odd>(0)을 사용하는 것을 알 수 있다. 이로써 수행 시간을 조금 줄일 수 있다. 또한 여기서 고속 푸리에 변환의 계산 효율성을 더 높일 수 있다. 위 그림에서 짝수와 홀수로 4개씩 이산 푸리에 변환을 수행했다. 나뉜 짝수, 홀수 이산 푸리에 변환들도 다시 그 안에서 짝수 홀수로 나뉠 수 있다. 예를 들어 위에서 4개의 짝수 이산 푸리에 변환은 홀수 번째에 있는 f(0)과 f(4)로 묶이고 짝수 번째에 있는 f(2)와 f(6)이 묶여서 그 안에서 다시 크기가 2인 짝수, 홀수 이산 푸리에 변환으로 만들 수 있다. 이렇게 반씩 계속 입력신호를 나눌 수 있는 데 이를 위해서 신호의 원소 수를 2의 거듭제곱 형태라고 가정한 것이다. 이렇게 재귀적으로 매 단계에서 변환을 분리하다 보면, 결국 최소 단위인 두 원소에 대한 변환을 수행하는 것까지 이루어진다. 재귀적인 분리를 마친 후에는 다음 그림과 같이 하면 된다. 다음과 같이 수행하면 각각 분리된 원소를 결합하는 과정이 logn이 되어 수행시간이 nlogn에 비례하게 된다. 또한 위 그림을 보면 결과 값을 구하는 과정을 도식한 선이 서로 교차한다고 해서 버터플라이 알고리즘이라고도 불린다.
위 그림을 보면 입력 신호의 순서가 순차적으로 되있지 않다는 것을 알 수 있다. 따라서 2의 거듭제곱 길이를 갖는 임의의 신호 배열이 위와 같은 순서가 되도록 만들어야 한느데 이를 스크램블이라고 한다. 스크램블은 원소의 원래 순서를 2진수로 표현하고 2진수 비트들의 순서를 거꾸로 한 후 이를 다시 10진수로 표현하여 수행하면 된다. 다음은 1차원 신호의 고속 푸리에 변환 메서드이다. 입력 신호의 배열 g[]가 주어졌을 때 g[1]이 첫 번째 원소의 실수부에 해당하고 g[2]가 첫 번째 원ㅅ의 허수부에 해당한다. isign은 변환이 정방향인지 역방향인지를 가리킨다.
void _FFT1d(double* g, unsigned long N, int isign)
{
int mmax, m, istep;
double wtemp, wr, wpr, wpi, wi, theta;
double tempr, tempi;
double temp;
// 스크램블 수행
int n = N*2;
int j = 1;
for (int i=1 ; i<n ; i+=2)
{
if (j > i)
{
temp = g[j];
g[j] = g[i];
g[i] = temp;
temp = g[j+1];
g[j+1] = g[i+1];
g[i+1] = temp;
}
m = n >> 1;
while (j>m && m>=2)
{
j -= m;
m >>= 1;
}
j += m;
}
// 버터플라이 알고리즘 수행
mmax = 2; // 2원소 DFT로 시작
while (n > mmax)
{
istep = mmax << 1;
theta = isign*(6.28318530717959/mmax);
wtemp = sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi = sin(theta);
wr = 1.0;
wi = 0.0;
for (m=1 ; m<mmax ; m+=2)
{
for (int i=m ; i<=n ; i+=istep)
{
j = i + mmax;
tempr = double (wr*g[j]-wi*g[j+1]);
tempi = double (wr*g[j+1]+wi*g[j]);
g[j] = g[i]-tempr;
g[j+1] = g[i+1]-tempi;
g[i] += tempr;
g[i+1] += tempi;
}
wr = (wtemp=wr)*wpr-wi*wpi+wr;
wi = wi*wpr+wtemp*wpi+wi;
}
mmax = istep;
}
if(isign == -1)
{
for (int i=1 ; i<= n ; i++)
{
g[i] /= N;
}
}
}※ 2차원 영상에 대한 고속 푸리에 변환
2차원 영상에 대한 고속 푸리에 변환은 1차원 변환을 가로 방향으로 수행하고 결과에 대하여 다시 세로 방향으로 변환을 수행하면 된다. 우리가 사용하는 신호 데이터의 배열은 행 우선 배열이기 때문에 세로 방향으로 푸리에 변환을 수행하려면 영상 배열의 가로와 세로 방향을 바꾸어 전치 행렬을 구하는 작업을 수행해야 한다. 다음 코드는 전치 영상을 구하는 Transpose() 메서드이다. Transpose() 메서드는 결과 영상으로 입력 영상의 가로, 세로 크기가 뒤바뀐 영상을 할당하고 입력 영상의 픽셀을 행 원으로 접근하면서 읽어들여 이를 결과 영상에 열 우선으로 할당한다. 고속 푸리에 변환에서는 입력 영사이 실수부와 허수부 두 채널을 가지기 때문에 여러 채널을 가진 영상에 대해서도 동작하게 만들었다.
CMyImage<T> Transpose() const
{
// 반환할 결과 영상 생성
CMyImage<T> dst(m_nHeight, m_nWidth, m_nChannels);
T* pSrc = m_pImageData;
T* pDst = dst.GetPtr();
int nWStepS = GetWStep();
int nWStepD = dst.GetWStep();
int nBack = (nWStepS*m_nHeight - m_nChannels);
for (int r=0 ; r<m_nHeight ; r++)
{
int idx = 0;
for (int c=0 ; c<m_nWidth ; c++)
{
// 각 채널 값을 순서대로 대입
for (int h=0 ; h<m_nChannels ; h++)
{
pDst[h] = pSrc[idx+h];
}
idx += m_nChannels; // 원본 영상은 오른쪽 픽셀로 이동
pDst += nWStepD; // 결과 영상은 아래쪽 픽셀로 이동
}
pSrc += nWStepS; // 원본 영상은 다음 행으로 이동
pDst -= nBack; // 결과 영상은 다음 열로 이동
}
return dst;
}복소수로 이루어지는 입력과 결과 신호는 실수부와 허수부를 분리하거나 결합해야 한다. 따라서 앞 포스트 "영상의 채널 단위 접근과 변환"에서 배운 PutChannelImg()와 GetChannelImg() 메서드를 사용할 것이다. 다음은 2차원 영상에 대한 고속 푸리에 변환 메서드이다.
CDoubleImage FFT(const CDoubleImage& imageIn)
{
int nWidth = imageIn.GetWidth();
int nHeight = imageIn.GetHeight();
//영상을 2의 거듭제곱으로 채우기 위해
int nFFTX = GetFFTSize(nWidth);
int nFFTY = GetFFTSize(nHeight);
//결과 영상을 0으로 초기화하고 실수부 채널에 입력영상 대입
CDoubleImage dst(nFFTX, nFFTY, 2);
dst.SetConstValue(0);
dst.PutChannelImg(imageIn, 0);
//행마다 1차원 FFT 수행
for (int r=0 ; r <nFFTY ; r++)
{
double* pDst = dst.GetPtr(r);
_FFT1d(pDst - 1, nFFTX, 1);
}
dst = dst.Transpose();
for (int r=0 ; r <nFFTX ; r++)
{
double* pDst = dst.GetPtr(r);
_FFT1d(pDst - 1, nFFTY, 1);
}
return dst.Transpose();
}
int GetFFTSize(int n)
{
return (int)pow(2.0, int(log((double)(n))/log(2.0) + 0.9999));
}
GetFFTSize 함수에서 log((doule)(n))/log(2.0)을 하게 되면 입력 받은 n이 2의 몇 승인지 알 수 있다. 여기에 1과 아주 가까운 수를 더해서 2의 pow 연산을 해주면 n이 넘는 2의 거듭제곱 가운데 가장 작은 수가 반환된다.
결과 영상은 복소수와 실수부로 구성되기 때문에 채널 수는 2로 잡는다. 또한 영상이 8btye double 형에 채널 수가 2개이기 때문에 행당 바이트 수가 4의 배수로 영상의 가로 크기와 같다.
다음은 역방향 고속 푸리에 변환 메서드를 볼 것이다. 주파수 분석 영상은 가로 세로 크기가 푸리에 변환을 통해 2의 거듭 제곱 형태로 변하기 때문에 역방향 고속 푸리에 변환 메서드에서는 가로와 세로 크기를 2의 거듭 제곱으로 맞추려고 호출하는 GetFFTSize 메서드를 호출하지 않아도 된다.
CDoubleImage IFFT(const CDoubleImage& imageIn)
{
int nFFTX = imageIn.GetWidth();
int nFFTY = imageIn.GetHeight();
CDoubleImage dst = imageIn;
for (int r=0 ; r <nFFTY ; r++)
{
double* pDst = dst.GetPtr(r);
_FFT1d(pDst - 1, nFFTX, -1);
}
dst = dst.Transpose();
for (int r=0 ; r <nFFTX ; r++)
{
double* pDst = dst.GetPtr(r);
_FFT1d(pDst - 1, nFFTY, -1);
}
return dst.Transpose();
}다음 사진은 FFT와 IFFT를 수행한 결과 영상이다.
'영상처리 프로그래밍' 카테고리의 다른 글
| 주파수 영역 영상 필터/가우시안 주파수 필터 (1) | 2017.06.22 |
|---|---|
| 주파수 스펙트럼과 위상각의 시각화 (1) | 2017.06.22 |
| SIFT 원리, 구현 (0) | 2017.06.19 |
| NCC, 교차 검토, 특징 기술자 (0) | 2017.06.19 |
| 영상 정합 (0) | 2017.06.19 |