#8 OpenCL 튜토리얼. Median Filter와 Bitonic Sort

Median FIlter

 Median Filter는 픽셀 주변을 둘러싼 픽셀 중 중간값을 찾아서 대체해주는 필터 입니다. 극대값과 극소값을 제거할 수 있습니다. 신호처리 나 이미지 처리에 자주 쓰이는 것으로 대표적으로는 이미지의 소금&후추 잡음을 없애는데 많이 쓰입니다. 다음 예에서는 이미지의 극대값을 없애는데 사용 하였습니다. 

Median Filter 연산

  Median FIlter연산은 연산의 대상이 되는 픽셀 주변의 픽셀을 한줄로 세워 가운데 서있는 픽셀의 값을 취하는 형태입니다. 다시말해 정렬을 해야합니다. 다들 아시다시피 정렬 연산의 Optimal은 O(n*logn)입니다. 수만 픽셀에 대해 정렬을 수행하는 것이라면 GPU를 이용하면 Optimal은 O(logn)입니다. 하지만 Median FIlter연산의 경우 각 픽셀에 해당하는 스레드, 다시말해 단일 프로세서에서 연산을 해야 합니다. 따라서 아무리 빠르게 연산을해도 O(n*logn)이며, 주변의 픽셀이라고 해봐야 커널의 크기가 13x13인경우를 생각해봐도 200개를 넘지 않습니다. 따라서 필요한 알고리즘은 단일 프로세서에서 비교적 적은 수의 항목에 대한 빠른 정렬을 가능하게 해주는 알고리즘 입니다. 

Bitonic sort

  Bitonic 정렬 알고리즘은 증가하는 시퀸스의 항목과 감소하는 시퀸스의 항목을 재귀적으로 교환하면서 결과적으로 증가하는 형태의 시퀸스를 만들 수있는 알고리즘입니다. Radix sort(기수정렬)이나 Merge sort(합병정렬)과 같은 경우에도 빠른 시간에 정렬이 가능한 알고리즘이지만 기수정렬의 경우에는 이미지 픽셀의 범위가 정수일때 적용 할 수 있기때문에 주로 Float값을 사용하는 이미지 처리에서는 사용하기 쉽지 않습니다. 합병정렬의 경우 GPU에서 사용하는 최적화 기법인 Reduction을 기법으로 응용하여 사용이 가능하지만 기본적으로 O(n*logn)의 복잡도를 가지는 정렬 알고리즘입니다. 

 Bitonic sort의 경우 (log^2(n))의 복잡도를 가집니다. 

Bitonic Sort의 연산 방법

   8개의 항목에 대해 Bitonic Sort하는 방법을 알아봄으로써 Bitonic Sort에 대해 알아보도록 하겠습니다. 4개의 줄에 대해 각각 설명 하도록 하겠습니다. 

1. 먼저 두개의 항목을 묶어 각각 오름차순, 내림차순으로 정렬합니다. 

2. 4개의 항목을 묶어 홀수항은 홀수항끼리 오름차순 짝수항은 짝수항끼리 오름차순. 반대편은 홀수항은 홀수항끼리 내림차순 짝수항은 짝수항끼리 내림차순 순으로 정렬합니다. 

3. 2개의 항목끼리 다시 오름차순, 내림차순이 되도록 정렬합니다. 

4. 최종적으로 오름차순인 시퀸스와 내림차순인 시퀸스가 완성됩니다. 

이후 왼쪽에서 4개의 덩어리에서 가장 큰 항과 오른쪽 4개의 덩어리에서 가장 작은 항을 비교하여 정렬해나가면 정렬이 완료됩니다. logn깊이만큼 연산하며 각 깊이마다 logn번의 연산을 합니다. 그리고, 이 연산은 GPU에서 아주 간편하게 표현이 가능합니다. 

  

  만약 항목의 크기가 짝수가 아니라면 어떻게 할까요? 만약 9개라면 위 방법을 적용하기전에 한쪽에 최대값을 비교하여 몰아넣고 이후연산에서 제외하면 됩니다. 

Bitonic Sort GPU 코드

  3x3 미디언 필터를 사용하는 예제를 이용해서 설명하도록 하겠습니다. GPU에서 Bitonic Sort를 표현하는 것은 아주 간단합니다. 위의 과정을 그대로 코드로 옮기는것이 끝입니다. 재귀적으로 적용되는 알고리즘이지만 최적화를 위해 한줄 한줄 옮겼기떄문에 매우 깁니다. 하지만 간단합니다. 

14~54 : 인덱스를 정하고 Global Memory에서 픽셀의 상하좌우, 대각선방향의점들의 위치를 계산하여 값을 uiRGBA[9]배열에 집어넣습니다. 이후 r0~ r8값에 값을 집어넣습니다.

56~135 : 두 값을 비교하며 최대값 최소값을 비교하여 위치를 정해줍니다. 위 그림에서 본 오름차순, 내림차순으로 정렬하는 과정입니다. GPU에서의 Min, Max연산은 비트 비교 연산으로 매우 빠르게 진행됩니다. 

136~147 : 최종적인 비교를통해 정렬을 마무리하고 중간값을 뽑습니다. 

//미디언필터
__kernel void MedianFIlter(__global float *inputImage, 
                            __global float*outputImage, 
                            int rows, 
                            int cols,
                             __global uchar *map, 
                             int orgWidth,
                             int orgHeight,
                              int startCol, 
                              int startRow)
{
 

    uint globalRow = get_global_id(1) ;
    uint globalCol = get_global_id(0) ;
    int dstIndex = (globalRow+startRow) * cols + globalCol+startCol;
 
    if(map[globalRow * orgWidth + globalCol] == 1){ 
    const int x = globalCol + startCol;
    const int y = globalRow + startRow;
 
    const int width = cols;
    const int iOffset = y * width;
    const int iPrev = iOffset - width;
    const int iNext = iOffset + width;
 
    uint uiRGBA[9];
 
    // get pixels within aperture
    uiRGBA[0] = inputImage[iPrev + x - 1];
    uiRGBA[1] = inputImage[iPrev + x];
    uiRGBA[2] = inputImage[iPrev + x + 1];
 
    uiRGBA[3] = inputImage[iOffset + x - 1];
    uiRGBA[4] = inputImage[iOffset + x];
    uiRGBA[5] = inputImage[iOffset + x + 1];
 
    uiRGBA[6] = inputImage[iNext + x - 1];
    uiRGBA[7] = inputImage[iNext + x];
    uiRGBA[8] = inputImage[iNext + x + 1];
 
    uint uiResult = 0;
    
    // extract next color channel
    uint r0,r1,r2,r3,r4,r5,r6,r7,r8;
    r0=uiRGBA[0];
    r1=uiRGBA[1];
    r2=uiRGBA[2];
    r3=uiRGBA[3];
    r4=uiRGBA[4];
    r5=uiRGBA[5];
    r6=uiRGBA[6];
    r7=uiRGBA[7];
    r8=uiRGBA[8];
 
    // perform partial bitonic sort to find current channel median
    uint uiMin = min(r0, r1);
    uint uiMax = max(r0, r1);
    r0 = uiMin;
    r1 = uiMax;
 
    uiMin = min(r3, r2);
    uiMax = max(r3, r2);
    r3 = uiMin;
    r2 = uiMax;
 
    uiMin = min(r2, r0);
    uiMax = max(r2, r0);
    r2 = uiMin;
    r0 = uiMax;
 
    uiMin = min(r3, r1);
    uiMax = max(r3, r1);
    r3 = uiMin;
    r1 = uiMax;
 
    uiMin = min(r1, r0);
    uiMax = max(r1, r0);
    r1 = uiMin;
    r0 = uiMax;
 
    uiMin = min(r3, r2);
    uiMax = max(r3, r2);
    r3 = uiMin;
    r2 = uiMax;
 
    uiMin = min(r5, r4);
    uiMax = max(r5, r4);
    r5 = uiMin;
    r4 = uiMax;
 
    uiMin = min(r7, r8);
    uiMax = max(r7, r8);
    r7 = uiMin;
    r8 = uiMax;
 
    uiMin = min(r6, r8);
    uiMax = max(r6, r8);
    r6 = uiMin;
    r8 = uiMax;
 
    uiMin = min(r6, r7);
    uiMax = max(r6, r7);
    r6 = uiMin;
    r7 = uiMax;
 
    uiMin = min(r4, r8);
    uiMax = max(r4, r8);
    r4 = uiMin;
    r8 = uiMax;
 
    uiMin = min(r4, r6);
    uiMax = max(r4, r6);
    r4 = uiMin;
    r6 = uiMax;
 
    uiMin = min(r5, r7);
    uiMax = max(r5, r7);
    r5 = uiMin;
    r7 = uiMax;
 
    uiMin = min(r4, r5);
    uiMax = max(r4, r5);
    r4 = uiMin;
    r5 = uiMax;
 
    uiMin = min(r6, r7);
    uiMax = max(r6, r7);
    r6 = uiMin;
    r7 = uiMax;
 
    uiMin = min(r0, r8);
    uiMax = max(r0, r8);
    r0 = uiMin;
    r8 = uiMax;
 
    r4 = max(r0, r4);
    r5 = max(r1, r5);
 
    r6 = max(r2, r6);
    r7 = max(r3, r7);
 
    r4 = min(r4, r6);
    r5 = min(r5, r7);
 
    // store found median into result
    uiResult |= min(r4, r5);

  상당히 길어보이는 코드입니다. O(n*logn) 알고리즘을 사용하면 약 30줄 정도면 표현이 가능할꺼 같은데 이렇게 복잡해 보이는 코드를 사용해야 할까 의문이 드실 수 있습니다. 제 컴퓨터에서는 50개 이하의 항목에 대해 최대 성능을 발휘한다는 Insertion Sort에 비해 약 8배의 성능 향상을 보였습니다. 실제로 퀵소트나 기타 정렬들의 모든 비교연산을 코드로 표현하면 위 연산보다 적은 양의 코드로 표현할 수 없다는걸 알 수 있습니다. 약간 어셈블리 코드와 같은 느낌이 드는 코드지요, 위 코드는 컴파일했을때의 코드와 거의 차이가 없는 코드라고 보시면 됩니다.