QGIS 파이썬 콘솔을 이용하여 오츠(Otsu) 임계치 기법 구현하기

REMOTE SENSING

QGIS 파이썬 콘솔을 이용하여 오츠(Otsu) 임계치 기법 구현하기

유병혁 2017. 7. 16. 04:50

이번 글에서는 QGIS 파이썬 콘솔을 이용하여 오츠(Otsu) 임계치 기법을 구현해 보겠습니다.

오츠 임계치는 영상을 분할하는 임계치를 자동으로 결정하는 비매개변수 무감독분류기법입니다.

영상이 있을 때 배경과 객체를 이진 분류할 때 주로 사용(MABLAB의 graythresh 함수)되는데요,

일부 연구에서는 드론 RGB 식생지수 영상을 식생과 비식생 영역으로 분류하는데 사용되기도 합니다.

오츠 임계치에 관한 내용은 아래 원문을 참조하시기 바랍니다.

논문: Otsu, N., 1979. A threshold selection method from gray-level histogram. IEEE Transactions on Systems, Man, and Cybernetics 9, 62–66.

PDF: http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=4310076

먼저, 오츠 임계치를 적용할 래스터 레이어를 추가해 보겠습니다.

저는 우리나라 지역을 촬영한 PROBA-V 1km NDVI 위성영상을 선택하였습니다.

# 소스 파일과 레이어 이름 정의
srcFile = "D:/PROBA-V_1km/PROBA-V_1km_20160101.tif" # source file
lyrName = "20160101" # layer name
# 레이어 생성
from qgis.core import QgsRasterLayer
rasLayer = QgsRasterLayer(srcFile, lyrName)

선택한 래스터의 밴드 개수, 셀 크기, 밴드 폭과 높이를 확인해 볼까요?!

참고로, 이 래스터는 WGS 84(EPSG:4326) 좌표계로 정의되어 있습니다.

# 밴드 개수
rasLayer.bandCount()
1
# 셀 크기
rasLayer.rasterUnitsPerPixelX(), rasLayer.rasterUnitsPerPixelY()
(0.008928571400000003, 0.008928571399999996)
# 밴드 폭, 높이
rasLayer.width(), rasLayer.height()
(897, 1121)

NDVI 래스터 색상표를 적용하고, 레이어 추가합니다. 해당 코드는 아래 글에 정리되어 있습니다.

QGIS 파이썬 콘솔을 이용하여 래스터 색상표(color ramp) 적용하기 | http://blog.daum.net/geoscience/1106

이번에는 래스터 레이어를 넘파이의 배열 클래스로 반환하도록 하겠습니다.

PyQGIS에서는 GDAL을 NumPy에 연결하는 gdalnumeric 라이브러리를 제공합니다.

# 넘파이 배열로 반환
import gdalnumeric
srcArr = gdalnumeric.LoadFile(srcFile)
srcArr
array([[ 0.16      ,  0.156     ,  0.152     , ..., -0.06      ,
        -0.06      , -0.064     ],
       [ 0.156     ,  0.152     ,  0.14399999, ..., -0.06      ,
        -0.064     , -0.036     ],
       [ 0.156     ,  0.152     ,  0.148     , ..., -0.064     ,
        -0.048     , -0.036     ],
       ..., 
       [ 0.016     ,  0.004     ,  0.016     , ...,  0.068     ,
         0.036     ,  0.016     ],
       [ 0.016     ,  0.016     ,  0.012     , ...,  0.068     ,
         0.024     ,  0.        ],
       [ 0.012     ,  0.02      ,  0.016     , ...,  0.052     ,
         0.94      , -0.02      ]], dtype=float32)

이제 본 글의 주제인 오츠 임계치를 계산해보도록 하겠습니다. 오츠 임계치는 2개 클래스 내

분산을 최소화시키거나, 2개 클래스 간 분산을 최대화 시키는 기준으로 t값을 찾는 방법입니다.

# 오츠(Otsu) 임계치 계산
import numpy as np
cMax, t = 0, 0
for x in np.arange(-1, 1, 0.01):
    clsL = srcArr[np.where(srcArr < x)]
    clsH = srcArr[np.where(srcArr >= x)]
    wL = float(clsL.size) / srcArr.size
    wH = float(clsH.size) / srcArr.size
    meanL = np.mean(clsL)
    meanH = np.mean(clsH)
    cVal = wL * wH * (meanL - meanH) ** 2
    if (cVal > cMax):
    cMax, t = cVal, x
print t
0.25