[ 영상 처리 ] Part3-1. OpenCV Edge Detection(C++)

이번시간에는 Median Filter를 복습한 후 Kernel Convolution을 Sobel Filter로 확장시켜 적용해본다.

또한, cv::Canny()가 아닌 Canny Edge Detectio의 동작 과정을 직접 구현해본다.

아래 포스팅을 완벽하게 이해한 후 따라오면 비교적 쉽다.

https://jinwoo-jung.com/60

[ 영상 처리 ] Part2-3. OpenCV Mask Processing(C++)

 

---

 

Median Filter

Median Filtering은 지난시간에 진행한 것 처럼, Filter 내의 Pixel Intensity 중 Median Intensity를 찾아 가공하는 픽셀의 Intensity로 설정하는 비선형 Filter이다.

 

Median Filtering

 

기존의 구현에서는 Mat::at()을 적용하여 찾은 Median Intensity로 해당 픽셀의 Intensity를 설정하였지만, Mat::data를 활용하여 조금 다르게 구현 해 보자.

 

myMedian()

void myMedian(const Mat& src_img, Mat& dst_img, const Size& kn_size) {
    dst_img = Mat::zeros(src_img.size(), CV_8UC1);

    int wd = src_img.cols; // 너비
    int hg = src_img.rows; // 높이
    int kwd = kn_size.width; // 커널 너비
    int khg = kn_size.height; // 커널 높이
    int rad_w = kwd / 2; // 반경 너비
    int rad_h = khg / 2; // 반경 높이

    // Data를 Pointer로 가져옴
    uchar* src_data = (uchar*)src_img.data;
    uchar* dst_data = (uchar*)dst_img.data;

    // Median을 계산하기 위해 동적 할당
    // Mat::data처럼 Kernel 크기 만큼의 공간을 동적 할당
    float* table = new float[kwd * khg]; 
    float tmp;

    for (int c = rad_w; c < wd - rad_w; c++) {
        for (int r = rad_h; r < hg - rad_h; r++) {

            // Kernel 내의 Intensity를 저장할 임시 공간
            tmp = 0.0f;

            for (int kc = -rad_w; kc <= rad_w; kc++) {
                for (int kr = -rad_h; kr <= rad_h; kr++) {

                    // Kernel 내의 Pixel Intensity를 모두 Table에 저장함
                    // 이때, (c,r)이 가공하고자 하는 중심 픽셀이고, 접근하는 순서는 열별로 진행 0,0 -> 1,0 -> ... -> 1,0 -> ...
                    tmp = (float)src_data[(r + kr) * wd + (c + kc)];
                    
                    // Intensity를 table에 저장
                    table[(kr + rad_h) * kwd + (kc + rad_w)] = tmp;
                }
            }

            // Kernel 내의 Intensity의 Median을 찾는 과정 : 정렬 후 중심값
            sort(table, table + kwd * khg); 
            dst_data[r * wd + c] = (uchar)table[(kwd * khg) / 2];
        }
    }

    delete[] table; 
}

 

하나하나 살펴보자.

    int wd = src_img.cols; // 너비
    int hg = src_img.rows; // 높이
    int kwd = kn_size.width; // 커널 너비
    int khg = kn_size.height; // 커널 높이
    int rad_w = kwd / 2; // 반경 너비
    int rad_h = khg / 2; // 반경 높이

    // Data를 Pointer로 가져옴
    uchar* src_data = (uchar*)src_img.data;
    uchar* dst_data = (uchar*)dst_img.data;

    // Median을 계산하기 위해 동적 할당
    // Mat::data처럼 Kernel 크기 만큼의 공간을 동적 할당
    float* table = new float[kwd * khg]; 
    float tmp;

 

Mat::data를 사용하여 원본 이미지의 픽셀과 결과 이미지의 픽셀에 접근할 것이다. 또한, Kernel 내의 픽셀 Intensity를 저장하기 위한 배열 src_data를 동적할당 해 주었다. 이때, 그 크기는 $kwd * khg$로 커널의 크기만큼 생성하였다. 

 

따라서 table은 1x9 형태의 Vector 처럼 존재할 것이고, 접근하는 방법은 Mat::data를 통해 접근하는 것 처럼 $[ (Kernel Width) * (접근하고자 하는 열) + 접근하고자 하는 행]$으로 접근하면 된다.

 

 

    for (int c = rad_w; c < wd - rad_w; c++) {
        for (int r = rad_h; r < hg - rad_h; r++) {

            // Kernel 내의 Intensity를 저장할 임시 공간
            tmp = 0.0f;

            for (int kc = -rad_w; kc <= rad_w; kc++) {
                for (int kr = -rad_h; kr <= rad_h; kr++) {

                    // Kernel 내의 Pixel Intensity를 모두 Table에 저장함
                    // 이때, (c,r)이 가공하고자 하는 중심 픽셀이고, 접근하는 순서는 열별로 진행 0,0 -> 1,0 -> ... -> 1,0 -> ...
                    tmp = (float)src_data[(r + kr) * wd + (c + kc)];
                    
                    // Intensity를 table에 저장
                    table[(kr + rad_h) * kwd + (kc + rad_w)] = tmp;
                }
            }

            // Kernel 내의 Intensity의 Median을 찾는 과정 : 정렬 후 중심값
            sort(table, table + kwd * khg); 
            dst_data[r * wd + c] = (uchar)table[(kwd * khg) / 2];
        }
    }

 

가장자리 문제를 해결하기 위해 Pixel을 접근하는 이중 for 문에서 Index가 rad_w 즉, Kernel Width의 절반 Index부터 wd - rad_w 즉, 전체 이미지 Width에서 rad_w를 뺀 Index까지 접근하는지는 이전 포스팅에서 자세히 설명하였으므로 생략하겠다. Height  또한 동일한 방법으로 접근하면 된다.

 

tmp는 임시 저장 공간으로, Kerenel 내부의 픽셀의 Intensity가 저장된다. 이를 해당 위치에 맞는 table 배열에 저장한다.

 

tmp = (float)src_data[(r + kr) * wd + (c + kc)];

 

src_data는 Mat::data이기 때문에, 위와 같이 접근하는 것을 이해하고 넘어가자.

 

// Kernel 내의 Intensity의 Median을 찾는 과정 : 정렬 후 중심값
sort(table, table + kwd * khg); 
dst_data[r * wd + c] = (uchar)table[(kwd * khg) / 2];

 

유효한 범위 내의 픽셀 Intensity가 table에 저장되어 있으므로 오름차순으로 정렬 후 Median Intensity를 현재 가공중인 픽셀 (c,r)의 Intensity로 설정함을 확인할 수 있다. 

 

Result

 

결과는 다음과 같다. 현재 가장자리 문제를 해결하기 위해서 해당 부분을 접근하지 않았기 때문에 테두리가 검은색(Intenstiy = 0)임을 확인할 수 있다. 또한, 이미지 내부의 Salt and Pepper Noise가 제거됨을 확인할 수 있다.

 

Sobel Edge Detection

들어가기 전에 Sobel Edge Detection은 아래 포스팅에서 자세히 다뤘으니 한번 공부하고 진행하길 바란다.

https://jinwoo-jung.tistory.com/62

[ 영상 처리 ] Ch4. Edge Detection

 

 

Sobel Filter는 Correlation Filter 중 하나로, Filter와 유사한 상관관계를 가지는 이미지에 대해서 강한 Response를 보이는 Filter이다. 

 

 

Sobel Filter는 위와 같이 구현되어 있다. Edge를 찾기 위해서는 Gradient를 구해야 하고 이는 Discrete한 Image의 경우 차분으로 대신한다고 공부하였다. 따라서 가로 혹은 세로를 기준으로 차분을 계산할 수 있도록 Filter가 구현되어 있다. 

 

Sobel Filter를 이용하여 Convolution하면 위 결과와 같이 x,y 각 방향으로의 Edge를 추출할 수 있다. 이때, 각 방향에 대하여 Pixel의 Intensity 분포가 Filter의 반대 방향 즉, 왼쪽에서 오른쪽으로 커지는 Filter지만, 픽셀은 오른쪽에서 왼쪽으로 Intensity가 커지는 경우에도 높은 Response를 보이나 부호가 반대이다. 따라서 모든 Edge를 고려하기 위해 Convolution 과정에서 절댓값을 씌울 것이다.

 

직접 구현해보자.

 

myKernelConv()

void myKernelConv(const Mat& src_img, Mat& dst_img, const Mat& kn) {
    dst_img = Mat::zeros(src_img.size(), CV_8UC1);

    int wd = src_img.cols; int hg = src_img.rows;
    int kwd = kn.cols; int khg = kn.rows;
    int rad_w = kwd / 2; int rad_h = khg / 2;

    float* kn_data = (float*)kn.data;
    uchar* src_data = (uchar*)src_img.data;
    uchar* dst_data = (uchar*)dst_img.data;

    float wei, tmp, sum;

    // < 커널 연산시 (가장자리 제외) >
    for (int c = rad_w + 1; c < wd - rad_w; c++) {
        for (int r = rad_h + 1; r < hg - rad_h; r++) {
            tmp = 0.f;
            sum = 0.f;
            // < 커널 연산시 >
            for (int kc = -rad_w; kc <= rad_w; kc++) {
                for (int kr = -rad_h; kr <= rad_h; kr++) {
                    wei = (float)kn_data[(kr + rad_h) * kwd + (kc + rad_w)];    // 해당 위치의 Kernel Data를 가져옴
                    tmp += wei * (float)src_data[(r + kr) * wd + (c + kc)];     // Pixel Intensity와의 연산결과   
                    sum += wei;                                                 // Kernel Weight의 합 
                }
            }

            // 절댓값을 통해 양방향의 Edge 모두 추출
            if (sum != 0.f) tmp = abs(tmp) / sum; // 정규화 및 overflow 방지
            else tmp = abs(tmp);

            if (tmp > 255.f) tmp = 255.f; // overflow 방지

            // dst_data는 CV_8UC1 Intensity Resolution을 가짐
            dst_data[r * wd + c] = (uchar)tmp;
        }
    }
}

 

하나하나 살펴보자.

    dst_img = Mat::zeros(src_img.size(), CV_8UC1);

    int wd = src_img.cols; int hg = src_img.rows;
    int kwd = kn.cols; int khg = kn.rows;
    int rad_w = kwd / 2; int rad_h = khg / 2;

    float* kn_data = (float*)kn.data;
    uchar* src_data = (uchar*)src_img.data;
    uchar* dst_data = (uchar*)dst_img.data;

    float wei, tmp, sum;

 

인자로 받아오는 kernel kn 역시 Mat 객체이기 때문에 Mat::data를 통해 Kernel Weight에 접근할 예정이다.

 

// < 커널 연산시 (가장자리 제외) >
for (int c = rad_w; c < wd - rad_w; c++) {
    for (int r = rad_h; r < hg - rad_h; r++) {
        tmp = 0.f;
        sum = 0.f;
        // < 커널 연산시 >
        for (int kc = -rad_w; kc <= rad_w; kc++) {
            for (int kr = -rad_h; kr <= rad_h; kr++) {
                wei = (float)kn_data[(kr + rad_h) * kwd + (kc + rad_w)];    // 해당 위치의 Kernel Data를 가져옴
                tmp += wei * (float)src_data[(r + kr) * wd + (c + kc)];     // Pixel Intensity와의 연산결과   
                sum += wei;                                                 // Kernel Weight의 합 
            }
        }

        // 절댓값을 통해 양방향의 Edge 모두 추출
        if (sum != 0.f) tmp = abs(tmp) / sum; // 정규화 및 overflow 방지
        else tmp = abs(tmp);

        if (tmp > 255.f) tmp = 255.f; // overflow 방지

        // dst_data는 CV_8UC1 Intensity Resolution을 가짐
        dst_data[r * wd + c] = (uchar)tmp;
    }
}

 

기존에 구현한 MyConv3x3()과 동일한 구조이다. 단순히 Weight가 Float Datatype인 것 제외하고는 다른 것이 없다.

Sobel Filter 기반으로 Edge를 검출할 때, 픽셀의 Intensity가 Filter Weight에 의해 음수 혹은 255를 넘어가는 경우가 존재한다. 따라서 이를 조건문을 통해 방지 해 주고, Weight의 합이 1이 아닌 경우 Normalization을 진행하였다.

 

Sobel Filter를 적용하기 전, Gaussian을 5x5 Kernel, $\sigma = 100$으로 Filtering 한 뒤, Sobel Filtering을 진행한 결과 다음과 같다.

 

Sobel Filtering Result

 

적용한 Filter를 유추해 보면 횡방향의 Edge가 검출되었으므로 횡방향으로 차분을 구하는 Sobel Filter임을 유추할 수 있다. 실제 적용한 Filter는 다음과 같다.

 

결국 이미지 내에서 횡방향으로 Pixel의 Intensity가 급격히 변하는 Edge가 크게 Response하여 흰색으로 결과 이미지에 나타남을 확인할 수 있다.

 

Canny Edge Detection

Canny Edge Detection은 아래와 같이 여러 Process를 단계적으로 수행하면서 Edge를 검출한다.

특히, 4번 과정에서 Min, Max Threshold를 설정 함으로써 검출할 Edge를 조절할 수 있다. 이를 직접 구현해보자.

 

각 과정별로 하나하나 살펴보자

1. Noise Reduction

Mat src_img = imread("rock.png", 0);
if (!src_img.data) printf("No image data \n");

Mat dst_img1, dst_img2, dst_img3, dst_img4, dst_img5;

Mat blur_img;
GaussianBlur(src_img, blur_img, Size(3, 3), 1.5);

 

cv::GaussianBlur를 활용하였고, 이때 Kernel Size 3x3, $\sigma = 1.5$이다.

 

2. Intensity Gradient Calculation

Mat magX = Mat(src_img.rows, src_img.cols, CV_32F);
Mat magY = Mat(src_img.rows, src_img.cols, CV_32F);
Sobel(blur_img, magX, CV_32F, 1, 0, 3);
Sobel(blur_img, magY, CV_32F, 0, 1, 3);

Mat sum = Mat(src_img.rows, src_img.cols, CV_64F);
Mat prodX = Mat(src_img.rows, src_img.cols, CV_64F);
Mat prodY = Mat(src_img.rows, src_img.cols, CV_64F);
multiply(magX, magX, prodX);
multiply(magY, magY, prodY);
sum = prodX + prodY;
sqrt(sum, sum);

 

먼저 cv::Sobel을 활용해 x,y 각 방향으로의 Gradient를 계산한다. 이후, 우리가 실제 사용하는 Gradient Maginutude는 다음과 같다.

 

현재 magX, magY는 차분으로 계산된 $\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}$만 존재한다. 따라서 cv::multiply()를 통해 제곱한 결과를 prodX에 저장한다.

 

cv::multiply

#include <iostream>
#include <vector>
#include <algorithm>
#include <time.h>
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"

using namespace cv;
using namespace std;


int main()
{
	// Sobel Filter
	float kn_data[] = { 1.f, 0.f, -1.f,
	                1.f, 0.f, -1.f,
	                1.f, 0.f, -1.f };
	float kn_data2[] = { 2.f, 0.f, 1.f,
					2.f, 0.f, 1.f,
					2.f, 0.f, 1.f };
	
	Mat kn(Size(3, 3), CV_32FC1, kn_data);
	Mat kn2(Size(3, 3), CV_32FC1, kn_data2);

	Mat Result;

	multiply(kn, kn2, Result);

	float* data = (float*)Result.data;

	for (int i = 0; i < 3; i++)
	{
		cout << data[i] << " ";
	}
	cout << endl;
	for (int i = 0; i < 3; i++)
	{
		cout << data[i] << " ";
	}
	cout << endl;
	for (int i = 0; i < 3; i++)
	{
		cout << data[i] << " ";
	}

	return 0;
}

 

multiply의 결과는 다음과 같다.

 

결국 첫번째 인자와 두번째 인자로 받은 Mat 객체를 동일한 위치의 Value끼리 곱한 결과를 세번째 인자에 저장하는 함수이다. 따라서 우리는 각각의 제곱한 결과를 sum에 저장하고 sqrt()를 통해 루트를 씌워 줌으로써 Gradient Magnitude를 계산할 수 있다.

 

Original Image                                                 ->                                          Gradient Magnitude

 

 

Mat magnitude = sum.clone();

Mat slopes = Mat(src_img.rows, src_img.cols, CV_32F);
divide(magY, magX, slopes);

 

또한, Gradient의 방향을 설정하기 위해서 cv::divide를 통해 나눠준다. 

 

우리는 Gradient의 크기와 방향을 모두 계산하였다. 이제는 Local Maximum을 찾기 위해 NMS를 적용한다.

nonMaximumSuppression()

void nonMaximumSuppression(Mat& magnitudeImage, Mat& directionImage) {
    Mat checkImage = Mat(magnitudeImage.rows, magnitudeImage.cols, CV_8U);
    MatIterator_<float >itMag = magnitudeImage.begin<float>();
    MatIterator_<float >itDirection = directionImage.begin<float>();
    MatIterator_<unsigned char>itRet = checkImage.begin<unsigned char >();
    MatIterator_<float > itEnd = magnitudeImage.end<float>();

    for (; itMag != itEnd; ++itDirection, ++itRet, ++itMag) {
        const Point pos = itRet.pos();
        float currentDirection = atan(*itDirection) * (180 / 3.142);
        while (currentDirection < 0) currentDirection += 180;
        *itDirection = currentDirection;
        if (currentDirection > 22.5 && currentDirection <= 67.5) {
            if (pos.y > 0 && pos.x > 0 && *itMag <= magnitudeImage.at<float>(pos.y - 1, pos.x - 1)) {
                magnitudeImage.at<float >(pos.y, pos.x) = 0;
            }
            if (pos.y < magnitudeImage.rows - 1 && pos.x < magnitudeImage.cols - 1 && *itMag <= magnitudeImage.at<float >(pos.y + 1, pos.x + 1)) {
                magnitudeImage.at<float>(pos.y, pos.x) = 0;
            }
        }
        else if (currentDirection > 67.5 && currentDirection <= 112.5) {
            if (pos.y > 0 && *itMag <= magnitudeImage.at<float >(pos.y - 1, pos.x)) {
                magnitudeImage.at<float>(pos.y, pos.x) = 0;
            }
            if (pos.y < magnitudeImage.rows - 1 && *itMag <= magnitudeImage.at<float >(pos.y + 1, pos.x)) {
                magnitudeImage.at<float>(pos.y, pos.x) = 0;
            }
        }
        else if (currentDirection > 112.5 && currentDirection <= 157.5) {
            if (pos.y > 0 && pos.x < magnitudeImage.cols - 1 && *itMag <= magnitudeImage.at<float >(pos.y - 1, pos.x + 1)) {
                magnitudeImage.at<float >(pos.y, pos.x) = 0;
            }
            if (pos.y < magnitudeImage.rows - 1 && pos.x>0 && *itMag <= magnitudeImage.at<float >(pos.y + 1, pos.x - 1)) {
                magnitudeImage.at<float >(pos.y, pos.x) = 0;
            }
        }
        else {
            if (pos.x > 0 && *itMag <= magnitudeImage.at<float >(pos.y, pos.x - 1)) {
                magnitudeImage.at<float >(pos.y, pos.x) = 0;
            }
            if (pos.x < magnitudeImage.cols - 1 && *itMag <= magnitudeImage.at<float >(pos.y, pos.x + 1)) {
                magnitudeImage.at<float >(pos.y, pos.x) = 0;
            }
        }
    }

}

 

조금 어려울 수 있지만, 하나하나 따라와 보자.

 

Mat checkImage = Mat(magnitudeImage.rows, magnitudeImage.cols, CV_8U);
MatIterator_<float >itMag = magnitudeImage.begin<float>();
MatIterator_<float >itDirection = directionImage.begin<float>();
MatIterator_<unsigned char>itRet = checkImage.begin<unsigned char >();
MatIterator_<float > itEnd = magnitudeImage.end<float>();

 

이전까지의 픽셀 접근 방법은 Matd::at()과 Mat::data를 사용하였다. 또 하나의 방법으로는 iterator에 접근하는 것이다.

---

잠깐 vector iterator에 대해 알아보자.

 

#include <iostream>
#include <vector>

using namespace std;

int main(void)
{
	vector<int> v;

	for (int i = 0; i < 7; i++) {
		v.push_back(10 * i);
	}

	vector<int>::iterator iter;

	iter = v.begin();

	cout << "Original" << endl;
	cout << &(*iter) << endl;
	cout << *iter << endl;

	iter += 2; // += 연산 사용
	cout << "After +2" << endl;
	cout << &(*iter) << endl;
	cout << *iter << endl;

	// 반복
	cout << "Every Element of vector v" << endl;
	for (iter = v.begin(); iter != v.end(); iter++) {
		cout << *iter << endl;
	}

	return 0;
}

 

iterator(반복자)는 Pointer 처럼 주소로 접근할 수 있다. 

iter = v.begin();

 

iterator iter에 vector의 시작 주소를 넣었다. 그러면 iter는 지금부터 vector의 시작 주소를 가르키는 pointer 이다. 따라서 해당 vector의 각 원소에 대한 접근은 다양한 방법으로 가능하다.

 

iter += 2;
*iter

iter[2];

 

pointer이기 때문에 2를 더하면 현재 vector의 Datatype은 int 형이기 때문에 한칸이 4byte를 차지하므로 8칸 뒤의 주소를 가르킨다. 따라서 Operator*로 해당 주소의 값을 나타낼 수 있다. 또한, 배열의 이름을 통해 접근하는 것처럼 $iter[2]$로도 동일하게 접근 가능하다. 

 

for (iter = v.begin(); iter != v.end(); iter++) {
    cout << *iter << endl;
}

 

iterator를 사용하는 가장 큰 이유는 다음과 같이 모든 배열에 접근할 때 iter에 배열의 시작 주소를 입력하고, 배열의 끝 주소 전까지 하나씩 증가시켜 가면서 모든 배열의 원소에 편리하게 접근할 수 있기 때문이다. 

---

Mat checkImage = Mat(magnitudeImage.rows, magnitudeImage.cols, CV_8U);
MatIterator_<float >itMag = magnitudeImage.begin<float>();
MatIterator_<float >itDirection = directionImage.begin<float>();
MatIterator_<unsigned char>itRet = checkImage.begin<unsigned char >();
MatIterator_<float > itEnd = magnitudeImage.end<float>();

 

다시 위 코드를 살펴보면 이해할 수 있을 것이다. 각 iterator는 magnitudeImage, directionImage, checkImage의 시작, 끝 주소를 가르킨다. 따라서 앞으로 iterator를 증가시켜 가면서 모든 값에 접근할 것이다. 

 

for (; itMag != itEnd; ++itDirection, ++itRet, ++itMag) 
{
    const Point pos = itRet.pos();
    float currentDirection = atan(*itDirection) * (180 / 3.142);
    while (currentDirection < 0) currentDirection += 180;
    *itDirection = currentDirection;
}

 

for문의 구조와 for문 안의 변수들에 대해서 먼저 살펴보자.

itMag != itEnd까지 반복 즉, magnitudeImage의 처음부터 끝까지 반복하는 것이다. 반복될 때 마다 magnitudeImage, directionImage, checkImage Iterator가 1씩 증가하여 다음 Pixel을 가르킴을 확인할 수 있다.

 

 

pos는 itRet.pos()이다. MaxIterator의 pose()는 현재 가르키고 있는 픽셀의 위치를 반환한다. 반환되는 Data Type은 Point로 멤버변수로 x,y를 가지며, 해당 멤버변수에 현재 itRet가 가르키는 위치가 반환된다.

 

currentDirection은 현재 Gradient의 Direction이다. Input으로 들어오는 directionImage에는 위에서 계산한 Gradient Direction이 각 픽셀의 Value로 저장되어 있다.

Mat slopes = Mat(src_img.rows, src_img.cols, CV_32F);
divide(magY, magX, slopes);

 

하지만 실제로 가르키는 방향의 각도를 구하기 위해서는 divide하여 저장된 값을 탄젠트 역함수인 atan()을 활용하여 각도로 변환하여야 한다. 또한, 이를 Radian으로 변환하기 위해 180을 곱하고 PI로 나눔을 확인할 수 있다.

 

 

한점에 대하여 여러 각도로 표현됨을 확인할 수 있는데, atan()은 절대각을 -90~90의 범위로 표현한다. 이때, Canny에서는 Direction을 결정하는데 0~180의 범위를 쓰기 때문에 만약 atan()으로 구한 각도가 0보다 작은 경우 180도를 더해서 해당 범위 내로 만들어준다. 이는 엣지 방향이 반대인 경우에도 하나의 방향으로 생각하기 때문이다. 이렇게 계산한 방향으로 원래 방향 값을 Update 한다.

 

각 방향에 대하여 NMS를 적용할 때, Canny의 경우 Direction을 크게 4 부분으로 나눈다. 이때, 픽셀 좌표계는 월드 좌표계와 달리 오른쪽 처럼 형성되었기 때문에 22.5~67.5의 범위가 아래와 같이 형성됨을 주의하자.

 if (currentDirection > 22.5 && currentDirection <= 67.5) {
     if (pos.y > 0 && pos.x > 0 && *itMag <= magnitudeImage.at<float>(pos.y - 1, pos.x - 1)) {
         magnitudeImage.at<float >(pos.y, pos.x) = 0;
     }
     if (pos.y < magnitudeImage.rows - 1 && pos.x < magnitudeImage.cols - 1 && *itMag <= magnitudeImage.at<float >(pos.y + 1, pos.x + 1)) {
         magnitudeImage.at<float>(pos.y, pos.x) = 0;
     }
 }

 

Gradient의 Direction이 22.5~67.5인 경우,  pos.y, pos.x가 0보다 큰 유효한 범위 내에서 해당 픽셀의 Gradient Magnitude가 대각선 위쪽 픽셀$(pos.x -1, pos.y -1)$의 Magnitude보다 작을 경우 해당 픽셀의 Magnitude를 0으로 설정한다. 또한, pos.y, po.x가 각각 Height -1, Width -1보다 작은 유효한 범위 내에서 해당 픽셀의 Gradient Magnitude가 대각선 아래쪽 픽셀$(pos.x +1, pos.y +1)$의 Magnitude보다 작은 경우 해당 픽셀의 Magnitude를 0으로 설정한다.

 

 

위 그림을 보면 조금 더 이해하기 쉽다. 4방향으로 나눈 구역에 해당하는 방향으로 Gradient Direction이 설정 되었으면, 해당 방향으로의 인접 픽셀에 대하여 최대값을 제외하고 모두 0으로 설정한다. 이 작업이 모든 픽셀에 대해서 모든 방향으로 진행되기 때문에 최종적으로 Local Max Gradient Magnitude만 남고 나머지는 0이 된다. 따라서 Edge의 두께를 줄일 수 있다.

 

nonMaximumSuppression()

void nonMaximumSuppression(Mat& magnitudeImage, Mat& directionImage) {
    Mat checkImage = Mat(magnitudeImage.rows, magnitudeImage.cols, CV_8U);
    MatIterator_<float >itMag = magnitudeImage.begin<float>();
    MatIterator_<float >itDirection = directionImage.begin<float>();
    MatIterator_<unsigned char>itRet = checkImage.begin<unsigned char >();
    MatIterator_<float > itEnd = magnitudeImage.end<float>();

    for (; itMag != itEnd; ++itDirection, ++itRet, ++itMag) {
        const Point pos = itRet.pos();
        float currentDirection = atan(*itDirection) * (180 / 3.142);
        while (currentDirection < 0) currentDirection += 180;
        *itDirection = currentDirection;
        if (currentDirection > 22.5 && currentDirection <= 67.5) {
            if (pos.y > 0 && pos.x > 0 && *itMag <= magnitudeImage.at<float>(pos.y - 1, pos.x - 1)) {
                magnitudeImage.at<float >(pos.y, pos.x) = 0;
            }
            if (pos.y < magnitudeImage.rows - 1 && pos.x < magnitudeImage.cols - 1 && *itMag <= magnitudeImage.at<float >(pos.y + 1, pos.x + 1)) {
                magnitudeImage.at<float>(pos.y, pos.x) = 0;
            }
        }
        else if (currentDirection > 67.5 && currentDirection <= 112.5) {
            if (pos.y > 0 && *itMag <= magnitudeImage.at<float >(pos.y - 1, pos.x)) {
                magnitudeImage.at<float>(pos.y, pos.x) = 0;
            }
            if (pos.y < magnitudeImage.rows - 1 && *itMag <= magnitudeImage.at<float >(pos.y + 1, pos.x)) {
                magnitudeImage.at<float>(pos.y, pos.x) = 0;
            }
        }
        else if (currentDirection > 112.5 && currentDirection <= 157.5) {
            if (pos.y > 0 && pos.x < magnitudeImage.cols - 1 && *itMag <= magnitudeImage.at<float >(pos.y - 1, pos.x + 1)) {
                magnitudeImage.at<float >(pos.y, pos.x) = 0;
            }
            if (pos.y < magnitudeImage.rows - 1 && pos.x>0 && *itMag <= magnitudeImage.at<float >(pos.y + 1, pos.x - 1)) {
                magnitudeImage.at<float >(pos.y, pos.x) = 0;
            }
        }
        else {
            if (pos.x > 0 && *itMag <= magnitudeImage.at<float >(pos.y, pos.x - 1)) {
                magnitudeImage.at<float >(pos.y, pos.x) = 0;
            }
            if (pos.x < magnitudeImage.cols - 1 && *itMag <= magnitudeImage.at<float >(pos.y, pos.x + 1)) {
                magnitudeImage.at<float >(pos.y, pos.x) = 0;
            }
        }
    }

}

 

4부분으로 나눈 모든 방향에 대해서 똫같은 매커니즘으로 NMS를 수행한다. 그러면 magnitudeImage로 받아온 Magnitude 가 저장된 원본에는 최댓값만 제외하고 모두 0으로 설정되어 있을 것이다.

 

Befor NMS                                                  ->                                                  After NMS

 

이제는 Edge를 추출해야 한다.

 

4. Double Thresholding
edgeDetect

void edgeDetect(Mat& magnitude, int tUpper, int tLower, Mat& edges) {
    int rows = magnitude.rows;
    int cols = magnitude.cols;

    edges = Mat(magnitude.size(), CV_32F, 0.0);

    // 픽셀 에지따기
    for (int x = 0; x < cols; x++) {
        for (int y = 0; y < rows; y++) {
            if (magnitude.at<float>(y, x) > tUpper) {
                followEdges(x, y, magnitude, tUpper, tLower, edges);
            }
        }
    }
}

 

Canny Edge Detection에서는 Min, Max Threshold이 존재한다. $T_{min}$보다 Gradient Magnitude가 작으면 Edge라고 판단하지 않고, $T_{min}$ 보다 크면 Edge라고 판단한다. $T_{min} ~ T_{max}$사이의 Magnitude가 존재하면 인접한 Magnitude를 따라가서 최종적으로 $T_{min}$보다 크게 되면 Edge라고 판단하게 된다.

 

따라서 코드를 보면 tUpper가 $T_{max}$이 되는데, 모든 픽셀에 대하여 픽셀의 Magnitude가 $T_{max}$보다 크면은 followEdges()를 동작시킨다. 이는 $T_{min} ~ T_{max}$사이의 Magnitude를 갖는 Edge들을 찾기 위함이다.

 

followEdges()

void followEdges(int x, int y, Mat& magnitude, int tUpper, int tLower, Mat& edges) {
    edges.at<float>(y, x) = 255;

    for (int i = -1; i <= 1; i++) {
        for (int j = -1; j <= 1; j++) {
            if ((i == 0) && (j == 0) && (x + i >= 0) && (y + j >= 0) &&
                (x + i < magnitude.cols) && (y + j < magnitude.rows)) {
                if ((magnitude.at<float>(y + j, x + i) > tLower) &&
                    (edges.at<float>(y + j, x + i) != 255)) {
                    followEdges(x + i, y + j, magnitude, tUpper, tLower, edges);
                }
            }
        }
    }
}

 

먼저 followEdges()의 인자로 온 픽셀 좌표의 Magnitude를 255로 만들어 Edge임을 명확히 한다.

 

이중 for문의 index를 보면 3x3 Kernel을 통해 인접픽셀에 접근한 것과 같이 현재 픽셀을 기준으로 각 방향의 인접 픽셀에 접근하는 것을 확인할 수 있다. 재귀적으로 함수가 호출되기 때문에 조금 복잡하지만 집중해보자.

 

 if ((i == 0) && (j == 0) && (x + i >= 0) && (y + j >= 0) && (x + i < magnitude.cols) && (y + j < magnitude.rows))

 

$ i == 0 , j == 0$ 이기 때문에 자기 자신을 의미한다. $ x+i >= 0 , y+j >= 0$ 즉, 유효한 범위 내에서 ..?

코드가 이상한데 내일 조교님께 문의 해보자.

 

즉, 인접 픽셀의 Magnitude가 $T_{min}$보다 크고, 255가 아닌 경우 즉, 이미 방문했던 픽셀이 아닌 경우 해당 픽셀도 Edge로 연결해야 하기 때문에 followEdges를 재귀적으로 호출해서 Magnitude를 255로 만들어 Edge로 판단할 수 있도록 한다. 

 

 

$T_{min} = 100, T_{max} = 200$으로 설정한 결과 오른쪽과 같이 Edge가 검출됨을 확인할 수 있다.

 

 

오른쪽 상단은 $T_{min} = 100, T_{max} = 200$, 왼쪽 하단은 $T_{min} = 100, T_{max} = 150$이다. $ T_{max}$가 작아질수록 Gradient Magnitude가 작은 Edge도 Edge라고 판단한다. 오른쪽 하단은 $T_{min} = 100, T_{max} = 50$ $인데, T_{max}$가 작아지니 이전엔 Edge가 아니라고 Fitering 된 Edge들도 검출되는 것을 확인할 수 있다.