집밖은 위험해

오늘 좌표계 때문에 손놓고 있던 허프 변환 문제를

겨우 겨우 해결해서 돌릴수가 있엇는데

내가 그동한 구현한 알고리즘들을

실제 사용하기에는 실행 속도가 너무 느렸다.

그러던 중에 이 알고리즘 설명, 구현 코드있던 사이트에서

픽셀 단위로 루프를 돌리는게 아닌

행렬 처리로 구하는 코드도 같이 구현해 두었더라.

그래서 한번 내가 보면서 구현한 코드랑

opencv 속도

벡터화 된 알고리즘 속도를 비교해보았다.

1. 픽셀 단위 처리를 통한 허프 라인 검출

에지 영상을 이미지로 주어 검출해보았는데

100 x 100 이미지임에도 2초 가까이 걸렷다.

opencv에선 어떨가

2. opencv houghline 함수 사용한 경우

위와 동일한 이미지에 opencv 함수만 사용해서 구해본 결과

200ms 정도 걸렸다.

위에서 구현한것보다 1/10 속도

3. 행렬 처리를 통한 허프 라인 검출

ref : towardsdatascience.com/lines-detection-with-hough-transform-84020b3b1549

이 사이트에서 소개하는 코드를 필요 없는 부분을 빼고 돌릴수 있게 고쳤다.

고친 결과 opencv와 동일한 속도를 보이더라

위 시험을 위한 수정 코드

def line_detection_vectorized(image, edge_image, num_rhos=180, num_thetas=180, t_count=220):
    res = image.copy()
    edge_height, edge_width = edge_image.shape[:2]
    edge_height_half, edge_width_half = edge_height / 2, edge_width / 2

    d = np.sqrt(np.square(edge_height) + np.square(edge_width))
    dtheta = 180 / num_thetas
    drho = (2 * d) / num_rhos

    thetas = np.arange(0, 180, step=dtheta)
    rhos = np.arange(-d, d, step=drho)

    cos_thetas = np.cos(np.deg2rad(thetas))
    sin_thetas = np.sin(np.deg2rad(thetas))

    accumulator = np.zeros((len(rhos), len(rhos)))

    edge_points = np.argwhere(edge_image != 0)
    edge_points = edge_points - np.array([[edge_height_half, edge_width_half]])
    #
    rho_values = np.matmul(edge_points, np.array([sin_thetas, cos_thetas]))

    accumulator, theta_vals, rho_vals = np.histogram2d(
      np.tile(thetas, rho_values.shape[0]),
      rho_values.ravel(),
      bins=[thetas, rhos]
    )
    accumulator = np.transpose(accumulator)
    lines = np.argwhere(accumulator > t_count)

    for line in lines:
        y, x = line
        rho = rhos[y]
        theta = thetas[x]
        a = np.cos(np.deg2rad(theta))
        b = np.sin(np.deg2rad(theta))
        x0 = (a * rho) + edge_width_half
        y0 = (b * rho) + edge_height_half
        x1 = int(x0 + 1000 * (-b))
        y1 = int(y0 + 1000 * (a))
        x2 = int(x0 - 1000 * (-b))
        y2 = int(y0 - 1000 * (a))
        res = cv2.line(res, (x1, y1), (x2, y2), (0, 255, 0), 1)

    plt.imshow(res,cmap="gray")
    return accumulator, rhos, thetas

일단 이 문제를 행렬 연산으로 풀수 있으면 

상당히 성능 개선할수 있는점을 볼수 있었다.

수정한 행렬 연산 코드를 나눠서 살펴보면

처음 대부분은 같은데

np.argwhere(edge_img != 0)을 사용한다.

np.argwhere(조건) 함수는 조건에 해당하는 인덱스를 반환해주더라

theta_resolution = 1
rho_resolution = 1
edge_img = morph.copy()
edge_height, edge_width = edge_img.shape
edge_height_half, edge_width_half = edge_height/2, edge_width/2
d = np.sqrt(np.square(edge_height) + np.square(edge_width))

thetas = np.arange(0, 180, theta_resolution)
rhos = np.arange(-d, d, rho_resolution)
cos_thetas = np.cos(np.deg2rad(thetas))
sin_thetas = np.sin(np.deg2rad(thetas))

accumulator = np.zeros((len(rhos), len(thetas)))
edge_points = np.argwhere(edge_img != 0)

edge_point는 조건에 맞는 인덱스 값들로 n x 2의 행렬 형태로 되어있다.

이미지가 100 x 100으로 1만 개의 픽셀이라면

edge_point는 706 x 2로 706개의 픽셀을 찾았다.

그다음에 

edge_points - 가로/2, 세로/2해주는데

이 연산으로 

이미지를 중심점(h/2, w/2)을 중심으로 위치했던 요소들이

원점(0,0)을 중심으로 이동하게 된다.

그 다음

edge_point와 [sin_thetas, cos_thetas]의 행렬 곱 연산을 수행한다.

이 의미를 천천히 보면

sin_thetas와 cos_thetas는 길이가 180인 벡터로

np.array()로 묶어주면

2 x 180의 행렬이 된다.

edge_points는

706 x 2형태의 행렬로 y, x 위치를 가지고 있으니

이 둘을 행렬 곱을 하면 706 x 180의 형태가 되며

행렬 곱의 결과는

존재하는 모든 에지 포인트의 y, x와 모든 범위의 s, c를 곱 한 후 더한 것으로 rho가 나오게 된다.

정리하면 모든 경우의 rho를 구하게 된다.

행이 i 번째 에지 포인트

열은 j 번째 theta_index

이 행렬[i, j] = rho 값이 된다.

그 다음에는 히스토그램2d 함수로

누적결과와 theta_vals, rho_vals를 받고 있는데

우선 내부 내용들 부터 살펴보자

np.tile은 처음 보는 함수인대

np.tile(a, b)

- a 배열을 b만큼 복사를 해준다고 한다.

ref : numpy.org/doc/stable/reference/generated/numpy.tile.html

thetas는 180, rho_values 708 x 180이므로

np.tile로 구한 값은 708 x 180 = 127080이 된다.

np.ravel()

- 주어진 배열을 1차원으로 평평 하게 만들어 준다.

ref : numpy.org/doc/stable/reference/generated/numpy.ravel.html

rho_values.ravel()로 모든 값들을 평활화 시켜주게 된다.

그다음 histogram2d를 보면

np.histogram2d

-numpy.histogram2d(x, y, bins=10, range=None, normed=None, weights=None, density=None)

- x : x vals

- y : y vals

- bins : binding vals

ref : numpy.org/doc/stable/reference/generated/numpy.histogram2d.html

이제 누적기를 위한 np.histogram2d를 수행한 결과

허프 공간이 나오고

5개의 점이 가장 진하게 보인다.

허프 공간을 구했으니

남은건 임계화와 시각화

정리 - 행렬 연산을 통한 허프 라인 검출

def hough_line(edge_img, rho_resolution= 1, theta_resolution = 1, threshold=10):
    edge_height, edge_width = edge_img.shape
    edge_height_half, edge_width_half = edge_height/2, edge_width/2
    d = np.sqrt(np.square(edge_height) + np.square(edge_width))
    
    thetas = np.arange(0, 180, theta_resolution)
    rhos = np.arange(-d, d, rho_resolution)
    cos_thetas = np.cos(np.deg2rad(thetas))
    sin_thetas = np.sin(np.deg2rad(thetas))
    
    accumulator = np.zeros((len(rhos), len(thetas)))
    edge_points = np.argwhere(edge_img != 0)
    
    edge_points = edge_points - np.array([[edge_height_half, edge_width_half]])
    rho_values = np.matmul(edge_points, np.array([sin_thetas, cos_thetas]))

    accumulator, theta_vals, rho_vals = np.histogram2d(
      np.tile(thetas, rho_values.shape[0]),
      rho_values.ravel(),
      bins=[thetas, rhos]
    )
    accumulator = np.transpose(accumulator)
    lines = np.argwhere(accumulator > threshold)

    res = img.copy()
    for line in lines:
        y, x = line
        rho = rhos[y]
        theta = thetas[x]
        a = np.cos(np.deg2rad(theta))
        b = np.sin(np.deg2rad(theta))
        x0 = (a * rho) + edge_width_half
        y0 = (b * rho) + edge_height_half
        x1 = int(x0 + 1000 * (-b))
        y1 = int(y0 + 1000 * (a))
        x2 = int(x0 - 1000 * (-b))
        y2 = int(y0 - 1000 * (a))
        res = cv2.line(res, (x1, y1), (x2, y2), (0, 255, 0), 1)

    plt.imshow(res,cmap="gray")
    return accumulator, rhos, thetas

ref : towardsdatascience.com/lines-detection-with-hough-transform-84020b3b1549

허프 변환을 이용한 직선 검출 알고리즘

0. 에지 영상 준비

1. rho와 theta의 범위 지정

- theta는 -90 ~ 90도

- rho는 -d ~ d. d는 이미지 대각선 길이.

- theta와 rho를 계산할 수 있도록 이산 값이 되게 해야함

2. 허프 공간을 나타낼 2차원 누적 배열 생성

- shape : n_rho * n_theta

- 0으로 초기화

- 가능한 모든 경우의 rho, theta들을 누적할 예정

3. 에지 영상의 모든 픽셀 확인. 에지 픽셀에 대해서 가능한 모든 theta와 rho를 구하여 허프 행렬에 해당 인덱스 증감

4. 지정한 임계치를 넘는 값의 인덱스(rho, theta)가 검출된 직선

알고리즘 구현 자체는 

여러 링크에 있는 내용들을 참고하다보니

대강 형태는 금방 만들었지만

시각화 할때가 많이 복잡했다.

이전 글에서 일반적인 xy 좌표계 상에서 rho, theta를 놓으면 아래 그림과 같이 그릴수 있지만

이미지 영역의 경우 y가 반대로 되어있는데

이때 rho, theta를 어떻게 보는 지가 문제였다.

제대로 이해하지 않은 상태에서 하다보니

rho와 theta를 xy 좌표계상에서 라인 플로팅을 할때,

rho와 theta로 그린 라인이 내가 생각했던 대로가 아니거나, 실제 에지와 일치 하지 않게 계속 나오더라

일단 사용한 기본 영상은 그림판에다가 일직선을 그려서 만들어봤다.

가우시안 블러링과 캐니 에지를 만들어주고

아무 생각없이(안해도 되지만) 모폴로지 닫힘도 해주었다.

이제 이 에지 영상에서 부터 시작하면 된다.

hough 공간의 누적 행렬을 만드는 함수를 구현하였다.

여기서 누적을 하는 이유는 많이 사용되는(임계치를 넘는) rho와 theta를 찾아 직선을 그려주기 위함인데,

하나의 에지 픽셀에 대해서 모든 방향 360도 -> 직선으로 나타낼태니 180도를 1도 간격(해상도)으로 나누어

각 에지의 180도 rho를 구하고

해당 rho값 자체(rho)와 theta의 인덱스(theta_idx)를

누적 행렬의 행과 열 인덱스로 사용하고 1씩 증감을 시켰다.

accumulator[rho, theta_idx] += 1

def hough_acc(edge_img, rho_resolution= 1, theta_resolution = 1):
    rows, cols = edge_img.shape    
    d = np.round(np.sqrt(rows**2 + cols**2))
    thetas = np.arange(-90, 90, theta_resolution)
    rhos = np.arange(-d, d + 1, rho_resolution)
    
    #all cos, sin thetas 0 from 180 degree
    cos_thetas = np.cos(np.deg2rad(thetas))
    sin_thetas = np.sin(np.deg2rad(thetas))
    accumulator = np.zeros((len(rhos), len(thetas)))
    cnt = 0
    for row in range(rows):
        for col in range(cols):
            # if seletecd point is edge
            if edge_img[row, col] != 0:
                cnt += 1
                # check all possible thetas
                for theta_idx in range(len(thetas)):
                    # rho = (x_dist * cos(theta)) + (y_dist * sin(theta))
                    rho =int(d + (col* cos_thetas[theta_idx])+(row *sin_thetas[theta_idx]))
                    accumulator[rho, theta_idx] += 1
    return accumulator, rhos, thetas    

위 함수로부터 acc와 rhos, theta를 전달받아 임계치를 넘는

rho와 theta들을 모아 반환해주는 함수를 구현하였다.

예를 들자면 아래와 같이 점 3개가 있는데

각 점에서 모든 방향으로 누적을 시키게 된다.

하지만 세점이 이어진 직선이 세 번 겹치므로 accumulator에서 가장 큰 값을 가지게 된다.

이 때, accumulator[rho, theta_idx] = 겹친 횟수이며

이 겹친 횟수가 임계치를 넘으면 직선으로 포함시킨다.

def find_peak(accumulator, rhos, thetas, threhsold=60):
    
    lines = []
    for y in range(accumulator.shape[0]):
        rho = rhos[y]
        for x in range(accumulator.shape[1]):
            if accumulator[y][x] > threhsold:
                theta = np.deg2rad(thetas[x])
                lines.append([rho, theta])
    
    return np.array(lines)

마지막으로 원본 이미지와

위 함수에서 반환한 직선으로 

라인을 플로팅 시키는 함수

rho와 theta로 부터

직선을 구하는데 

x0, y0가 (그리고자 하는) 법선과 만나는 점이고,

x1,y1 x2,y2는 이미지 밖에서 부터 그려주기 위해 구한 좌표가 된다.

def draw_hough_lines(img, lines):
    res = img.copy()
    for i, line in enumerate(lines):
        rho = line[0]
        theta = line[1]
        # reverse engineer lines from rhos and thetas
        a = np.cos(theta)
        b = np.sin(theta)
        x0 = a*rho
        y0 = b*rho
        
        # these are then scaled so that the lines go off the edges of the image
        x1 = int(x0 + 1000*(-b))
        y1 = int(y0 + 1000*(a))
        x2 = int(x0 - 1000*(-b))
        y2 = int(y0 - 1000*(a))
        cv2.circle(res, (int(x0), int(y0)),2,(0,0,255),1)
        cv2.line(res, (x1, y1), (x2, y2), (0, 255, 0), 1)
        plt.imshow(res,cmap="gray")

위 이미지로 부터 에지를 구할수가 있었다.

acc, rhos, thetas = hough_acc(morph)
lines = find_peak(acc, rhos,thetas,threhsold=40)
print(lines[:5])
draw_hough_lines(img, lines)

그런데 구한 rho와 theta 값이 이해되지 않았다.

rho = 38, theta= -0.7853인데

theta는 라디안이므로 육십분법으로 변환하면

- 44가 되더라

어딜 기준으로 -44도로 38만큼 떨어진 지점에서 직선을 그렷길래 찾았나 했더니

이미지 평면의 원점이 기준이더라.

일반 적인 xy평면에서

y 양의 방향으로 theta가 +, y 음의 방향으로 theta가 -가 됬었는데

y방향이 반대로 되면서 theta가 +인 경우 아래로, theta가 -인 경우 위로 가는거더라

이부분 때문에 좀 많이 햇갈렸다.

원래는 일반적인 사진에서 직선 검출을 하려고 했지만

픽셀 단위로 연산하고, 모든 픽셀들을 다루다보니 상당히 시간이 오래 걸리더라

이 부분을 보니 아에 행렬, 벡터 단위로 처리하는 코드도 있더라

그렇게 수정하면 조금 나아 지지않을까 싶다.

대신 위 이미지에 라인을 몇개 더 그려서 검출 해봤다.

임계치를 40으로 주었을때 아래와 같은 결과를 얻을수 있었다.

여기서 lines의 앞 다섯개 요소를 보니

rho가 -65, theta가 -1.38 정도 되더라

이게 어느직선인가 보기 위해서

호도법을 육십분법으로 바꾸고 원점에서 찾아보았다.

rho = 65, theta= - 80 정도로 대강 저 에지에서 검출한 직선인것 같다.

ref : towardsdatascience.com/lines-detection-with-hough-transform-84020b3b1549

허프 변환 Hough Transform

- 파울 V.C 허프가 1962년 사진에 존재하는 직선들을 찾아내기 위해 제안한 알고리즘.

- 허프 변환은 컴퓨터 비전 분야에서 직선이나 원 같은 특징을 검출할때 주로 사용.

- 에지 영상(주로 케니 에지)에서 허프 변환을 통해 허프 공간(a, b 공간 표현법: 아래에서 후술)상의 교점을 구함.

- 허프 공간의 교점은 xy 공간에서 직선으로 표현됨

기본 허프 변환

- 기존의 이미지는 x-y 평면 상에서 표현됨

- 위에서 말한 a, b 공간은 a(기울기), b(y 절편)로 이루어진 공간을 말함.

- x, y 공간에서 한 점은, ab 공간(허프 공간)에서 한 직선으로 표현할 수 있음.

- xy 공간 상에 점이 하나 더 생기면 a-b 공간상의 교점이 생긴다.

- a-b 공간 상의 교점이 두 점이 지나는 직선의 기울기와 y 절편 값이 된다.

 => 이렇게 하여 ab 표현을 이용한 직선 검출이 수행.

ab 공간의 문제

- 아래의 좌측과 같은 에지 영상이 주어지면, 이들은 대부분 ab 공간 상에서 하나의 교점에서 만날것임

- 하지만 에지가 수직인 경우, ab 공간에서 b가 무한대가 되버리는 문제 발생

 -> 대안으로 rho, phi를 이용한 법선 normal line 표현법을 사용함.

 * 법선 normal line : 곡선의 한 점을 지나며, 접선에 수직인 직선

법선을 이용한 허프 변환

- 두 점을 지나는 법선은 원점까지 수직인 거리가 rho, 이때 각을 phi로 표현

- 수직 거리 rho, 각도 phi의 관계 : rho = x cos(phi) + y sin(phi)

- 하단의 우측 그림처럼 이전과 마찬가지로 rho, phi 공간 상에서 교점이 xy 공간상에서 직선을 의미함.

- ab 공간떄와는 달리 에지가 수직 혹은 수평의 형태인 경우에도 rho나 phi가 무한대로 발산하지 않음.

정리한 내용은 얼마 안되는데 시간이 너무 늦어졌다...

구현은 내일로

이미지 그라디언트를 마치고 허프 변환을 어떻게 구현할까 생각해보고 있었는데,

opencv를보니 에지 영상에서 허프 변환을 수행했었다.

그런데 내가 만든 이미지 그라디언트는 커널 사이즈도 3 x 3으로 고정되어 있고,

그렇게 좋은 이미지들을 만드는것 같지는 않아보였다.

그래서

이미지 그라디언트 구현을 공부하면서 조금 연습했으니

이번에는 opencv 문서를 보면서 캐니 에지 검출기를 구현해 보려고 한다.

ref 

docs.opencv.org/master/da/d22/tutorial_py_canny.html

캐니 에지 검출기 canny edge detector

- 존 F 캐니가 만든 많이 쓰이는 에지 검출기

- 아래의 단계로 구성

캐니 에지 검출기 구성 단계

1. 노이즈 제거

- 에지 검출기는 노이즈에 강인해야 한다. 그러니 우선  5 x 5 가우시안 필터로 노이즈를 제거하자.

2. 이미지 그라디언트 강도, 방향 구하기

- 이미지를 가우시안 필터로 블러링 해주고 나서, 소벨 커널로 수평, 수직 방향으로 필터링해서 그라디언트 구하기

- 이미지 그라디언트를 연산해서 에지 그라디언트 강도와 방향을 구해주자.

* 그라디언트 방향은 에지에 수직.

* 에지 방향 각도는 아래의 식으로 구한 후 반올림을 하여 (0, 45, 90, 135도)  중 하나가 되도록 한다.

3. 비최대 억제 Non-maximu suppression

- 그라디언트 크기와 방향을 구했다면, 필요없는 픽셀들을 없앨 준비가 되었다.

- 그러려먼 우선 모든 픽셀들에 대해 에지 그라디언트 방향 부근의 픽셀들이 지역 최대치인지 확인하여야 한다.

- 아래 그림을 보면 점 A는 에지이고, 그라디언트 방향은 에지에 수직이며, 점 B와 C는 그라디언트 방향에 있다.

- 점 B와 C를 보고 A가 지역 최대점인지 판단한다. A가 지역 최대가 아니라면, 억제되어 0이 된다.

4. 이력 현상 임계화 hystersis thresholding(이력 현상 : 이전의 상태가 현재에 영향을 주는 현상)

- 이제 모든 에지들이 정말로 에지인지 아닌지 판단하기 위해서 두 임계값 minVal, maxVal이 필요하다.

- 그라디언트 강도가 max val보다 크면 에지이고, min val보다 작으면 에지가 아니라고 판단해서 버린다.

- 만약 임계치에 있는 경우에는 연결성을 고려해서 에지 여부를 판단한다.

- 아래의 에지 A는 maxVal보다 크니 확실히 에지로 보지만, 에지 C는 에지 A에 연결되었으니 에지로본다.

- 하지만 에지 B는 minVal보다 크고, 에지 C랑 같은 곳에 있지만 A랑 연결이 안되있으므로 에지가 아닌것으로 본다.

2020. 12. 1  1:50

- 위 이론 내용 정리 후 자료 찾아가면서 구현하다보니 늦어졌다.

- 코드 일부분, 시각화 결과를 우선 보인 후 전체 코드를 맨 뒤에다가 배치.

전체 구성

1. 기본 이미지 시각화

2. 가우시안 필터링 적용

3. 이미지 그라디언트 방향 영상, 그라디언트 크기 영상 시각화

4. 비최대 억제 결과 시각화

5. 이력 현상 임계화 추가 -> 캐니 에지 검출기 완성

6. 캐니 에지 검출 과정 모든 이미지 시각화

7. 전체 코드

후기

- 기존의 opencv에서 제공하는 함수에 비하면 원하는 결과도 잘 나오지 않고, 속도도 매우 느린 편이다.

- 하지만 직접 수식, 개념들을 완전히는 아니더라도 조금씩 이해해 가면서 구현하는 과정이 좋았다.

- 이걸 하는 중간에 보니 내가 연결성 부분이나 범람 채움 알고리즘 같은 내용들을 넘어갔던것 같다. 내일 한번 봐야할듯.

1. 기본 이미지 시각화

- opencv 예제에서 제공하는 기본 messi5.jpg

import numpy as np
import matplotlib.pyplot as plt
from utils import padding, gaussian_filtering
from PIL import Image
import cv2

plt.figure(figsize=(8,12))
img = Image.open("./res/messi5.jpg").convert("L")
img = np.asarray(img)
plt.imshow(img, cmap="gray")

2. 가우시안 필터링 적용

- opencv doc에서 kernel size가 5인 가우시안 커널을 사용한다고 한다. 시그마는 0.8로 지정

plt.figure(figsize=(8,12))
img_blur = gaussian_filtering(img,k_size=5, sigma=0.8)
plt.imshow(img_blur,cmap="gray")

3. 이미지 그라디언트 방향 영상, 그라디언트 크기 영상 시각화

grad_intensity = get_gradient_intensity(img_blur)
plt.title("gradient direction image")
plt.imshow(grad_intensity[:,:,0])

plt.title("gradient magnitude image")
plt.imshow(grad_intensity[:,:,1])

4. 비최대 억제 결과 시각화

- 그라디언드 방향에 존재하는 픽셀 세 개 비교.

- 중앙 픽셀이 가장 크면 나두고, 같은 방향 다른 픽셀이 크면 억제하여 0

supressed_img = non_maximum_suppression(grad_intensity)
plt.imshow(supressed_img)

5. 이력 현상 임계화 추가 -> 캐니 에지 검출기 완성

img = Image.open("./res/messi5.jpg").convert("L")
img = np.asarray(img)
img_blur = gaussian_filtering(img,k_size=5, sigma=1)
grad_intensity = get_gradient_intensity(img_blur)
suppressed_img = non_maximum_suppression(grad_intensity)
canny_edge = hysteresis_threshold(suppressed_img=suppressed_img,min_th=100, max_th=200)


plt.figure(figsize=(64,32))
plt.subplot(6, 1, 1)
plt.title("original")
plt.imshow(img)
plt.subplot(6, 1, 2)
plt.title("gaussian blur")
plt.imshow(img_blur)
plt.subplot(6, 1, 3)
plt.title("gradient direction")
plt.imshow(grad_intensity[:, :,0])
plt.subplot(6, 1, 4)
plt.title("gradient magnitude")
plt.imshow(grad_intensity[:, :,1])
plt.subplot(6, 1, 5)
plt.title("non maximum suppression ")
plt.imshow(suppressed_img)
plt.subplot(6, 1, 6)
plt.title("canny edge image")
plt.imshow(canny_edge)

6. 캐니 에지 검출 과정 모든 이미지 시각화

def canny_edge_visualize(img, min_th=100, max_th=200):
    """
    canny_edge_visualize
    - visualize all images from original to canny edge
    
    parameter
    - img : original image
    - min_th : minimal threshold value
    - max_th : maximum threshold value
    """
    img_blur = gaussian_filtering(img,k_size=5, sigma=1)
    grad_intensity = get_gradient_intensity(img_blur)
    suppressed_img = non_maximum_suppression(grad_intensity)
    canny_edge = hysteresis_threshold(suppressed_img=suppressed_img,min_th=100, max_th=200)

    plt.figure(figsize=(64,32))
    plt.subplot(6, 1, 1)
    plt.title("original")
    plt.imshow(img)
    plt.subplot(6, 1, 2)
    plt.title("gaussian blur")
    plt.imshow(img_blur)
    plt.subplot(6, 1, 3)
    plt.title("gradient magnitude")
    plt.imshow(grad_intensity[:, :,0])
    plt.subplot(6, 1, 4)
    plt.title("gradient direction")
    plt.imshow(grad_intensity[:, :,1])
    plt.subplot(6, 1, 5)
    plt.title("non maximum suppression ")
    plt.imshow(suppressed_img)
    plt.subplot(6, 1, 6)
    plt.title("canny edge image")
    plt.imshow(canny_edge)
   
img = Image.open("./res/kid.jpg").convert("L")
img = np.asarray(img)
canny_edge_visualize(img)

7. 전체 코드

"""
Canny Edge Detector
- 2020. 12.1 01:50

"""


"""
ref
- https://en.wikipedia.org/wiki/Canny_edge_detector
"""
def sobel_kerenl():
    kernel_x = np.array([
        [-1, 0, 1],
        [-2, 0, 2],
        [-1, 0, 1]
    ])
    kernel_y = np.array([
        [1, 2, 1],
        [0, 0, 0],
        [-1, -2, -1]
    ])
    return kernel_x, kernel_y

def get_rounded_gradient_angle(ang):
    """
    get_rounded_gradient_angle
    - gradient direction angle round to one of four angles
    
    return : one of four angles
            representing vertical, horizontal and two diagonal direction
    
    parameteres
    ---------------
    ang : gradient direction angle
    """
    vals = [0, np.pi*1/4, np.pi*2/4, np.pi*3/4, np.pi*4/4]
    interval = [np.pi*1/8, np.pi*3/8, np.pi * 5/8, np.pi* 7/8]

    if ang < interval[0] and ang >= interval[-1]:
        ang = vals[0]
    elif ang < interval[1] and ang >= interval[0]:
        ang = vals[1]
    elif ang < interval[2] and ang >= interval[1]:
        ang = vals[2]
    elif ang < interval[3] and ang >= interval[2]:
        ang = vals[3]
    else:
        ang = vals[4]
    return ang
        

def get_gradient_intensity(img):
    """
    get gradient_intensity
    - calculate gradient direction and magnitude
    
    return (rows, cols, 2) shape of image about grad direction and mag
    
    parameteres
    ------------
    img : blured image
    """
    k_size = 3
    rows, cols = img.shape
    kernel_x, kernel_y = sobel_kerenl()
    pad_img = padding(img, k_size=k_size)
    res = np.zeros((rows, cols, 2))
    sx, sy = 0, 0

    for i in range(0, rows):
        for j in range(0, cols):
            boundary = pad_img[i:i+k_size, j:j+k_size]
            sx = np.sum(kernel_x * boundary)
            sy = np.sum(kernel_y * boundary)            
            ang = np.arctan2(sy, sx)
            mag = abs(sx) + abs(sy)
            if ang < 0:
                ang = ang + np.pi
            ang = get_rounded_gradient_angle(ang)
            res[i, j, 0] = ang
            res[i, j, 1] = mag
    return res



def check_local_maximum(direction=None, boundary=None):
    """
    check_local_maximum
    - check if center value is local maximum depend on gradient direction
    return True if center valus is local maximum
    
    parameter
    ------------
    direction : gradient direction
    boundary : region of image for finding local maximum
    """
    if direction == 0: # 0 degree, east and west direction
        kernel = np.array([
            [0, 0, 0],
            [1, 1, 1],
            [0, 0, 0]
        ])

    elif direction == 1: #45 degree, north east and south west direction
        kernel = np.array([
            [0, 0, 1],
            [0, 1, 0],
            [1, 0, 0]
        ])
    elif direction == 2: #90 degree, north & south direction
        kernel = np.array([
            [0, 1, 0],
            [0, 1, 0],
            [0, 1, 0]
        ])
    else : #135 degree, north west & south east direction
        kernel = np.array([
            [1, 0, 0],
            [0, 1, 0],
            [0, 0, 1]
        ])
    
    max_val = np.max(kernel * boundary)
    if boundary[1,1] == max_val: #local maximum
        return True
    else: # not local maximu
        return False

    
    
def non_maximum_suppression(grad_intensity=grad_intensity):
    """
    non_maximum_suppression
    - check a pixel wiht it's neighbors in the direction of gradient
    - if a pixel is not local maximum, it is suppressed(means to be zero)
    """
    directions = [0, np.pi*1/4, np.pi*2/4, np.pi*3/4, np.pi*4/4]
    
    k_size = 3
    grad_direction = grad_intensity[:, :, 0]
    grad_magnitude = grad_intensity[:, :, 1]
    
    rows, cols = grad_magnitude.shape
    pad_img = padding(grad_magnitude, k_size=k_size)
    res_img = np.zeros((rows,cols))
    sx, sy = 0, 0
    
    
    for i in range(0, rows):
        for j in range(0, cols):
            direction = directions.index(grad_direction[i,j])
            boundary = pad_img[i:i+k_size, j:j+k_size]
            if check_local_maximum(direction, boundary) == True:
                res_img[i, j] = grad_magnitude[i, j]
            else:
                res_img[i, j] = 0
    return res_img



def hysteresis_threshold(suppressed_img=None, min_th=None, max_th=None):
    """
    hysterisis_threshold
    - check edge is connected to strong edge using eight connectivity
    
    parameter
    ---------------------
    - suppressed img : non maximum suppressed image
    - min_th : minimal threshold value
    - max_th : maximum threshold value
    """
    k_size = 3
    rows, cols = supressed_img.shape
    pad_img = padding(suppressed_img, k_size=k_size)
    res_img = np.zeros((rows,cols))

    eight_connectivity = np.array([
        [1, 1, 1],
        [1, 0, 1],
        [1, 1, 1]
    ])    
    
    for i in range(0, rows):
        for j in range(0, cols):
            if pad_img[i+1, j+1] < min_th:
                res_img[i, j] = 0
            else:
                boundary = pad_img[i:i+k_size, j:j+k_size]
                max_magnitude = np.max(boundary * eight_connectivity)
                if max_magnitude >= max_th : # pixel is connected to real edge
                    res_img[i, j] = 255
                else:
                    res_img[i, j] = 0
    return res_img


def canny_edge_detector(img=None, min_th=None, max_th=None):
    """
    canny_edge_detector
    - detect canny edge from original image
    
    parameter
    - img : original image
    - min_th : minimal threshold value
    - max_th : maximum threshold value
    """
    img_blur = gaussian_filtering(img,k_size=5, sigma=1)
    grad_intensity = get_gradient_intensity(img_blur)
    suppressed_img = non_maximum_suppression(grad_intensity)
    canny_edge = hysteresis_threshold(suppressed_img=suppressed_img,
                                      min_th=min_th, max_th=max_th)
    return canny_edge


def canny_edge_visualize(img, min_th=100, max_th=200):
    """
    canny_edge_visualize
    - visualize all images from original to canny edge
    
    parameter
    - min_th : minimal threshold value
    - max_th : maximum threshold value
    """
    img_blur = gaussian_filtering(img,k_size=5, sigma=1)
    grad_intensity = get_gradient_intensity(img_blur)
    suppressed_img = non_maximum_suppression(grad_intensity)
    canny_edge = hysteresis_threshold(suppressed_img=suppressed_img,min_th=100, max_th=200)

    plt.figure(figsize=(64,32))
    plt.subplot(6, 1, 1)
    plt.title("original")
    plt.imshow(img)
    plt.subplot(6, 1, 2)
    plt.title("gaussian blur")
    plt.imshow(img_blur)
    plt.subplot(6, 1, 3)
    plt.title("gradient magnitude")
    plt.imshow(grad_intensity[:, :,0])
    plt.subplot(6, 1, 4)
    plt.title("gradient direction")
    plt.imshow(grad_intensity[:, :,1])
    plt.subplot(6, 1, 5)
    plt.title("non maximum suppression ")
    plt.imshow(suppressed_img)
    plt.subplot(6, 1, 6)
    plt.title("canny edge image")
    plt.imshow(canny_edge)

    

이번에 뭘 다뤄야 하나 생각하다가

opencv document를 보니

이번에는 이미지 그라디언트를 볼 차래더라

이미지 그라디언트는 이미지의 기울기를 구하는 것을 말한다.

이 기울기는 영상내 급격히 변하는 곳으로 영상 내 물체의 윤곽선(에지)가 된다.

이에 대한 설명은 다음 링크와 cmu 슬라이드에서 잘 설명해주고 있다.

ref 

iskim3068.tistory.com/49

www.cs.cmu.edu/~16385/s17/Slides/4.0_Image_Gradients_and_Gradient_Filtering.pdf

그러면 이미지 그라디언트 코드를 구현해보자.

우선 소벨 필터를 이용한 이미지 그라디언트

"""
ref
- https://en.wikipedia.org/wiki/Edge_detection
- http://www.cs.cmu.edu/~16385/s17/Slides/4.0_Image_Gradients_and_Gradient_Filtering.pdf
- https://iskim3068.tistory.com/49
"""
def sobel_kerenl():
    kernel_x = np.array([
        [-1, 0, 1],
        [-2, 0, 2],
        [-1, 0, 1]
    ])
    kernel_y = np.array([
        [1, 2, 1],
        [0, 0, 0],
        [-1, -2, -1]
    ])
    return kernel_x, kernel_y

def sobel(img, method=None):
    """
    get image gradient using sobel operater
    
    parameteres
    ------------
    img : input image applying sobel filter
    method : 1(x direction), 2(y dicrection), 3(x + y direction)
    """
    k_size = 3
    rows, cols = img.shape
    kernel_x, kernel_y = sobel_kerenl()
    pad_img = padding(img, k_size=k_size)
    res_img = np.zeros((rows,cols))
    sx, sy = 0, 0

    for i in range(0, rows):
        for j in range(0, cols):
            boundary = pad_img[i:i+k_size, j:j+k_size]
            if method == 1:
                sx = np.sum(kernel_x * boundary)
            elif method == 2:
                sy = np.sum(kernel_y * boundary)
            else:
                sx = np.sum(kernel_x * boundary)
                sy = np.sum(kernel_y * boundary)
            res_img[i,j] = np.sqrt(sx**2 + sy**2)
    
    return res_img

x 방향, y 방향, x+y방향 영상 시각화 결과

plt.figure(figsize=(12,12))
for i in range(1, 4):
    plt.subplot(2, 2, i)
    sobel_img = sobel(img, method=i)
    plt.imshow(sobel_img,cmap="gray")

라플라시안 그라디언트 영상

def laplacian_filter():
    kernel_x = np.array([
        [0, 1, 0],
        [1, -4, 1],
        [0, 1, 0]
    ])
    kernel_y = np.array([
        [0, -1, 0],
        [-1, 4, -1],
        [0, -1, 0]
    ])
    return kernel_x, kernel_y


def laplacian(img):
    """
    get image gradient using laplacian filter
    
    parameteres
    ------------
    img : input image applying laplacian filter
    """
    k_size = 3
    rows, cols = img.shape
    kernel_x, kernel_y = laplacian_filter()
    pad_img = padding(img, k_size=k_size)
    res_img = np.zeros((rows,cols))
    sx, sy = 0, 0

    for i in range(0, rows):
        for j in range(0, cols):
            boundary = pad_img[i:i+k_size, j:j+k_size]
            sx = np.sum(kernel_x * boundary)
            sy = np.sum(kernel_y * boundary)
            res_img[i,j] = np.sqrt(sx**2 + sy**2)    
    return res_img
 
laplacian_img = laplacian(img)
plt.imshow(laplacian_img,cmap="gray")

모폴로지 연산

- 컴퓨터 비전의 영역 연산 중 한 종류

- 이미지를 깍거나 팽창시켜준다. -> 침식, 팽창 연산

- 침식 연산과 팽창 연산을 이용하여 열림 연산과 닫힘 연산이 가능.

- 열림 연산 : 침식 -> 팽창. 이미지를 깍아 노이즈를 제거, 깍여진 이미지를 팽창하여 원래 형태로 복원

- 닫힘 연산 : 팽창 ->침식. 이미지를 팽창시켜, 내부 구멍들을 채움, 팽창된 이미지를 침식하여 원래 형태로

1. 이미지 로드. 기본 이미지 시각화

import numpy as np
import matplotlib.pyplot as plt
from utils import Histogram, Threshold, padding
from PIL import Image


img = Image.open("./res/j.png").convert("L")
img = np.asarray(img)

plt.imshow(img, cmap="gray")

2. 모폴로지 연산 구현

def erosion(boundary=None, kernel=None):
    """
    erosion operation
    - a pixel element is 0 at least under kernel is 0
    -> return 0
    """
    boundary = boundary * kernel
    if (np.min(boundary) == 0):
        return 0
    else:
        return 255

def dilation(boundary=None, kernel=None):
    """
    erosion operation
    - a pixel element is not 0 at least under kernel is not 0
    -> return 255
    """
    boundary = boundary * kernel
    if (np.max(boundary) != 0):
        return 255
    else:
        return 0

def openning(img=None, k_size=None):
    """
    openning operation
    - erosion followed by dilation
    - it removes noise
    """
    erosion_img = morphology(img=img, method=1, k_size=k_size)
    opened_img = morphology(img=erosion_img, method=2, k_size=k_size)
    return opened_img

def closing(img=None, k_size=None):
    """
    closing operation
    - dilation follwed by erosion
    - it can close small holes inside the objects. 
    """
    dilation_img = morphology(img=img, method=2, k_size=k_size)
    closed_img = morphology(img=dilation_img, method=1, k_size=k_size)
    return closed_img


def morphology(img=None, method=None, k_size=None):
    """
    input
    img : input image
    method : 1(erosion), 2(dilation), 3(openning), 4(closing)
    k_size : kernel size
    
    output
    res_img : morphology operation image
    """
    rows, cols = img.shape
    pad_img = padding(img, k_size)
    kernel = np.ones((k_size,k_size))
    res_img = img.copy()
    
    if method == 1 or method == 2:
        for i in range(0, rows):
            for j in range(0, cols):
                if method == 1: #erosion operation
                    res_img[i, j] = erosion(pad_img[i:i+k_size, j:j+k_size], kernel)
                elif method == 2: #
                    res_img[i, j] = dilation(pad_img[i:i+k_size, j:j+k_size], kernel)
    if method == 3:
        res_img = openning(img, k_size=k_size)
    elif method == 4:
        res_img = closing(img, k_size=k_size)

    return res_img

3. 침식 연산 결과

res = morphology(img, method=1, k_size=5)
plt.title("erosion operation with k_size = 5")
plt.imshow(res, cmap="gray")

4. 팽창 연산 결과

res = morphology(img, method=2, k_size=5)
plt.title("dilation operation with k_size = 5")
plt.imshow(res, cmap="gray")

5. 열기 연산

5.1 열기 연산 이전. 노이즈 이미지

noised_j = Image.open("./res/noised_j.PNG").convert("L")
noised_j = np.asarray(noised_j)
plt.title("befor opening operation")
plt.imshow(noised_j, cmap="gray")

5.2 열기 연산 수행

res = morphology(noised_j, method=3, k_size=7)
plt.title("after opening operation with k_size :7")
plt.imshow(res, cmap="gray")

6. 닫기 연산

6.1 닫기 연산 전. 이미지

small_hole = Image.open("./res/j_small_hole.PNG").convert("L")
small_hole = np.asarray(small_hole)
plt.title("before closing operation")
plt.imshow(small_hole, cmap="gray")

6.2 닫기 연산 후

res = morphology(small_hole, method=4, k_size=5)
plt.title("after closing operation with k_size :7")
plt.imshow(res, cmap="gray")

스무딩을 왜 했었더라..

영상 잡음 노이즈 제거용으로 스무딩 처리를 한다고 하더라

가장 많이 사용되는게 가우시안 스무딩, 블러링이라고도 한다.

그런데 지금까지 내가 구현한 코드가 잘못된 부분이 있었다.

평균 스무딩, 평균 임계화를 다루면서

내가 패딩 처리를 제대로 하지 않은것도 있엇고

사실 커널, 마스크를 만들어 곱 연산 후 나누어주었어야 했는데, 그냥 평균으로 내버린 점이 그렇다.

일단 평균 임계/스무딩까지는 그렇게 해도 괜찬았는데

가우시안 스무딩을 구현하면서가 문제가 되더라..

ref : youngest-programming.tistory.com/235?category=877383

그러다 한 학생이 구현한 내용들을 봤는데

보기 좋고, 이해하기도 쉽고 정말 잘 만들었더라.

나는 대충 돌아가게만 만들어놓고 내버리거나 하는 경우가

아무튼 이분 내용 참고해서 간단하게 나마 구현할 수 있었다.

가우시안 커널, 마스크를 어떻게 구현하나 막막 했는데 많은 참고가 되었다.

def gaussian_kernel(k_size, sigma):
    """
    param
    k_size : Gaussian kernel size
    sigma : gaussian kernel standard variance
    
    return
    filter = k_size * k_size gaussian filter
    """
    size = k_size//2
    y, x = np.ogrid[-size:size+1, -size:size+1]
    #ref : https://en.wikipedia.org/wiki/Gaussian_filter
    filter = 1/(2*np.pi * (sigma**2)) * np.exp(-1 *(x**2 + y**2) /(2*(sigma**2)))
    sum = filter.sum()
    filter /= sum
    return filter

내가 이전에 만든 코드에는 인덱싱 최대, 최소 범위를 지정해서 하다보니

이런 마스크와 이미지 부분을 컨벌루션 시킬수가 없었다.

이 문제를 패딩 이미지를 만들어서 해결할수 있었다.

def padding(img, k_size):
    """
    param
    img : padding img
    k_size : kernel size
    
    return 
    res : padded img
    """
    pad_size = k_size//2
    rows, cols, ch = img.shape
    
    res = np.zeros((rows + (2*pad_size), cols+(2*pad_size), ch), dtype=np.float)
    
    if pad_size == 0:
        res = img.copy()
    else:
        res[pad_size:-pad_size, pad_size:-pad_size] = img.copy()
    return res
    

인덱싱은 기존 이미지 사이즈 대로 인덱싱대로 되지만

곱 연산은 패딩 이미지와 커널을 중심으로 수행된다.

커널을 기존 이미지에 적용하여 이미지 밖으로 넘어가는 문제는 제거됬다.

def gaussian_filtering(img, k_size=3,sigma=1):
    """
    param
    img : input img
    k_size : kernel size
    sigma : standard deviation
    
    return
    filtered_img : gaussian filtered image returned
    """
    rows, cols, channels = img.shape
    filter = gaussian_kernel(k_size, sigma)
    pad_img = padding(img,k_size)
    filtered_img = np.zeros((rows, cols, channels), dtype=np.float32)
    
    for ch in range(0, channels):
        for i in range(rows):
            for j in range(cols):
                filtered_img[i, j, ch] = np.sum(filter * pad_img[i:i+k_size, j:j+k_size, ch])

    return filtered_img.astype(np.uint8)

마지막으로 시각화

커널 사이즈가 작을때에는 sigma를 크게주어도 큰 차이가 없지만

사이즈가 클땐 sigma 값에 따라 블러링 정도가 강해지는걸 볼수 있다.

plt.figure(figsize=(12,12))
row = 3
col = 3

for i in range(0, 3):
    for j in range(0,3):
        plt.subplot(3, 3, 1+ i*3 + j)
        k_size = 3+2*i
        sigma = 3*j + 1
        res = gaussian_filtering(img, k_size = k_size, sigma = sigma)
        title = "k size : " + str(k_size) + ",  sigma = " + str(sigma)
        plt.title(title)
        plt.imshow(res)

컴퓨터 비전에서의 연산으로

점 연산, 영역 연산 등이 있었던걸로 기억한다.

스무딩은 아마 영역 연산이었던것 같은데

이전 코드를 조금 수정해서 평균 스무딩을 구현하였다.

먼저 기본 이미지 부터 읽어보고

컬러 이미지를 다룰수 있도록 기존의 Histogram을 조금 수정 후 플로팅 시켰다.

마지막은 커널 사이즈별 블러링 차이

import numpy as np
import matplotlib.pyplot as plt
from utils import Histogram, Threshold
from PIL import Image

img = Image.open("./res/lena.jpg").convert("RGB")
img = np.asarray(img)

plt.imshow(img)


def Histogram(img):
    row, col, channels = img.shape
    hist = np.zeros((256, channels))
    for channel in range(0,channels):
        for i in range(0, row):
            for j in range(0, col):
                hist[img[i, j, channel], channel] += 1
    
    return hist

hist = Histogram(img)
plt.figure(figsize=(12,12))
plt.subplot(3,1,1)
plt.bar(np.arange(0,256), hist[:,0])
plt.subplot(3,1,2)
plt.bar(np.arange(0,256), hist[:,1])
plt.subplot(3,1,3)
plt.bar(np.arange(0,256), hist[:,2])



def meanBlur(img,  block_size=5):

    if type(img) is not np.ndarray:
        raise AssertionError("img is not ndarray")
    row, col, channels = img.shape


    res = np.zeros((row, col, channels),dtype=int)
    if (block_size % 2 == 0):
        block_size += 1
    
    for ch in range(0, channels):
        for i in range(0, row):
            for j in range(0, col):
                x_min = j-block_size//2
                x_max = j+block_size//2
                y_min = i-block_size//2
                y_max = i+block_size//2

                if x_min <= 0:
                    x_min = 0

                if x_max >= col:
                    x_max = col

                if y_min <= 0:
                    y_min = 0

                if y_max >= row:
                    y_max = row


                val = img[y_min:y_max, x_min:x_max, ch].mean()
                res[i, j, ch] = int(val)
    return res
    
    
res = meanBlur(img, block_size=3)
plt.figure(figsize=(12,12))
plt.subplot(3,1,1)
res = meanBlur(img, block_size=3)
plt.imshow(res)
plt.subplot(3,1,2)
res = meanBlur(img, block_size=5)
plt.imshow(res)
plt.subplot(3,1,3)
res = meanBlur(img, block_size=9)
plt.imshow(res)

다음 링크에서 오츠 이진화에 대해 잘 정리하여 참고

medipixel.github.io/post/2019-06-07-image-thresholding/

* 본 내용을 정리하면서 앞의 일부 내용들은 엉망이라 유의하자.

오츠 이진화

- 최적의 임계치를 찾아 영상을 이진화하는 알고리즘

- 흑색 영상과 백색 영상의 분산 가중 합을 가장 적게하는 임계치가 최적 임계치

- 위 내용을 대강 그림으로 그리면 히스토그램에 대해서 다음과 같다

인줄 알았는데

해당 블로그에서 hat h(t)에 대한 설명부분이나, between class variance에 대한 설명이 빠져있어

조금 해매버렸다.

이 블로그도 조금은 잘못됬지만 대강

j07051.tistory.com/364

클래스내 분산 within class variance와 클래스간 분산 between class variance에 대해서 보다 잘 설명해주고 있다.

다만 틀린 부분은 

클래스 내 분산이 최소가 되도록 하는 임계치 t의 인덱스가 오츠 임계값이고,

이는 클래스간 분산이 최대가 되는것과 동일한 의미이다.

오츠 임계치를 다시 정리하면,

클래스 내 분산이 최소가 되는 인덱스나

클래스 간 분산이 최대가 되는 인덱스를 찾는 것이고

이 인덱스가 가장 잘 영상을 이진화 시키는 오츠 임계치 값이 된다.

이런 식으로 클래스 간 분산이 최대(클래스 내 분산이 최소)로 분할하는 임계치이다.

아래의 코드는 내가 구현한 오츠 임계치 함수인데

def otsuThreshold(img):
    
    if type(img) is not np.ndarray:
        raise AssertionError("img is not ndarray")
    
    hist = Histogram(img)
    
    vars_within = []
    vars_between = []

    zero = 1.e-17
    for t in range(0, 256):
        sumb = np.sum(hist[:t]) + zero
        sumw = np.sum(hist[t:]) + zero
        sum = sumb + sumw
        wb = sumb/sum
        ww = sumw/sum
        
        mub = zero
        muw = zero
        for i in range(0, t):
            mub += i * hist[i]/sumb
        for i in range(t, 256):
            muw += i * hist[i]/sumw

        vb = zero
        vw = zero
        for i in range(0, t):
            vb += hist[i] * ((i - mub)**2)/sumb
        for i in range(t, 256):
            vw += hist[i] * ((i - muw)**2)/sumw
    
        var_within = wb * vb + ww * vw
        vars_within.append(var_within)


    th = vars_within.index(min(vars_within))
    print(th)
    res = Threshold(img, th)
    return res

res = otsuThreshold(img)

위 함수를 사용한 결과 오츠 임계치는 96가 되며 아래의 이미지를 얻었다.

히스토그램으로 보면 가장 흑백 영상 클래스가 잘 분할하도록 분포하고 있다.

* 위 이미지 처리 관련 코드들은 사이킷런에서도 잘 구현되 있었다.

* 내것보다 유효성 검증, cumsum 같은 함수 사용, 인덱싱 처리에서 훨씬 간결하다.

ref : github.com/scikit-image/scikit-image/blob/master/skimage/filters/thresholding.py#L237

아무튼 이론과 구현에 대한 내용들은

learning opencv에서 잘 제공하고 있었다

www.learnopencv.com/otsu-thresholding-with-opencv/

적응 평균 임계치 이진화에 대한 내용은

opencv document와 다음 링크를 참고했다.

homepages.inf.ed.ac.uk/rbf/HIPR2/adpthrsh.htm

opencv-python.readthedocs.io/en/latest/doc/09.imageThresholding/imageThresholding.html

이전에 단순 이진화의 경우 본인이 지정한 밝기 값을 기준으로 이진화를 수행하였다.

평균 기반 적응적 임계치 이진화 같은 경우 본인이 임계치를 지정하는 것이 아니라

해당 픽셀을 중심으로 블록 사이즈 만큼의 커널(마스크) 영역의 밝기 평균을 지역 임계치로 사용한다.

해당 픽셀 밝기와 지역 임계치를 비교하여 255를 줄지 0을 줄지 판단한다.

다음 함수는 구현내용

def adaptiveThresholdMean(img,  block_size=5, C=4):

    if type(img) is not np.ndarray:
        raise AssertionError("img is not ndarray")
    row, col = img.shape


    res = np.zeros((row, col))
    if (block_size % 2 == 0):
        block_size += 1
    
    for i in range(0, row):
        for j in range(0, col):
            x_min = j-block_size//2
            x_max = j+block_size//2
            y_min = i-block_size//2
            y_max = i+block_size//2
            
            if x_min <= 0:
                x_min = 0
            
            if x_max >= col:
                x_max = col
            
            if y_min <= 0:
                y_min = 0
            
            if y_max >= row:
                y_max = row

            
            val = img[y_min:y_max, x_min:x_max].mean()
            local_th = val-C
            if img[i,j] >= local_th:
                res[i, j] = 255
            else:
                res[i, j] = 0
    return res

다음 예시는 아래의 링크에서 소개하는 설정대로 해보았다.

homepages.inf.ed.ac.uk/rbf/HIPR2/adpthrsh.htm

dave.jpg는 opencv 다큐먼트꺼