หาขอบที่เป็นเส้นตรงด้วย Probabilistic Hough Transform

นอกจาก Hough Transform แล้ว OpenCV ยังมีอีกคำสั่งหนึ่งในการหาขอบเส้นตรง นั่นคือ HoughLinesP ซึ่งสามารถศึกษาวิธีการได้จาก [1]

คำสั่ง HoughLinesP จะให้ผลลัพธ์เป็นจุดสองจุดที่เป็นต้น (x1,y1) และปลาย (x2,y2) ของเส้นตรงที่คำนวณได้ การพลอตผลลัพธ์ก็จะทำได้ง่ายขึ้น

ลองดูตัวอย่าง เปรียบเทียบกับบทความก่อนหน้าครับ

import cv2
import numpy as np

img = cv2.imread("building.jpg")
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
#Appy Gaussian blur to remove some noises
inoise = cv2.GaussianBlur(gray,(3,3),sigmaX=0)

#Canny
lowThresh = 50
upThresh = 200
edge = cv2.Canny(inoise,lowThresh,upThresh)

#Probabilistic Hough Line Transform
lines = cv2.HoughLinesP(edge,rho=1,theta=np.pi/180,threshold=120,minLineLength=30,maxLineGap=10)

img2 = img.copy()
#Plot detected lines
for x1,y1,x2,y2 in lines[0]:
    cv2.line(img2,(x1,y1),(x2,y2),(0,255,0),2)

cv2.imshow("Probabilistic Hough Lines",img2)
cv2.waitKey()
cv2.destroyAllWindows()

อ้างอิง
[1] Matas, J. and Galambos, C. and Kittler, J.V., Robust Detection of Lines Using the Progressive Probabilistic
Hough Transform. CVIU 78 1, pp 119-137 (2000)

การหาขอบที่เป็นเส้นตรงด้วย Hough Line Transform

เมื่อเราหาขอบ (edge) ของรูปภาพได้แล้ว เราจะทราบได้อย่างไรว่ามีขอบที่เป็นเส้นตรงหรือไม่

เราสามารถใช้ Hough Transform เพื่อตอบคำถามนี้ได้ ดังนี้

ถ้ามีเส้นตรงใดๆ (สีดำ) ตามรูป ที่มีระยะจากจุด (0,0) ไปตั้งฉากกับเส้นตรงนั้นเท่ากับระยะ r โดยเส้นตั้งฉากนี้ทำมุมกับแกน x เป็นมุม T

เส้นตรงปกติจะมีสมการเป็น y = mx + c

แต่เราสามารถเขียนอยู่ในรูปของ polar coordinate คือ
y = (-1/tan(T)) x + (r/sin(T))
y = -cos(T)/sin(T) x  + r/sin(T)
y sin(T) = -x cos(T) + r
หรือ
y sin(T) + x cos(T) = r

โดยทั่วไปแล้ว จะกำหนดค่า T ให้อยู่ในช่วง [0, PI] และค่า r ใน [-D, D] เมื่อ D คือความยาวของเส้นทะแยงมุมของภาพ (ระยะไกลสุดของภาพที่วัดจากจุดกำเนิด)

สมมติว่าเรามีจุดๆหนึ่ง เส้นตรงที่ลากผ่านจุดนี้จะมีได้ไม่จำกัด นั่นคือ จะมีค่า r และ T มากมาย

แต่ถ้ามีจุดสองจุดขึ้นไป เส้นตรงที่ลากผ่านจุดเหล่านี้จะมีได้เส้นเดียว นั่นคือ เส้นที่มีค่า r และ T เท่ากัน

ดังนั้น เมื่อเราคิดใน polar coordinate หลักการหาขอบที่เป็นเส้นตรงด้วย Hough Transform ก็คือ

1. หาขอบของภาพด้วยวิธีิใดๆ เช่น Canny ให้ได้ผลลัพธ์เป็นขาวดำ (ขาว = ขอบ)
2. เริ่มต้นด้วย r=0, T=0 (เส้นตรงที่ทับกับแกน y)  ไล่ตรวจสอบว่ามีพิกเซลสีขาว (ขอบ) ที่อยู่บนแนวเส้นตรงนี้เท่าใด ถ้ามีเกินกว่าค่า threshold ที่กำหนด (เช่น 100 จุด) ก็ให้เก็บค่า (r,T) นี้ไว้
3. เปลี่ยนแปลงค่า r ไปเรื่อยๆ ตรวจสอบเช่นเดิม
4. เปลี่ยนแปลงค่า T ไปเรื่อยๆ ตรวจสอบเช่นเดิม
5. สุดท้ายจะได้ชุดของ (r,T) ซึ่งแต่ละคู่จะเป็นค่าของขอบเส้นตรงที่อยู่ในรูป
6. พลอตเส้นตรงเหล่านั้น

ใน OpenCV เราสามารถใช้คำสั่ง cv2.HoughLines() ดังตัวอย่างต่อไปนี้

import cv2
import numpy as np

img = cv2.imread("building.jpg")
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
#Appy Gaussian blur to remove some noises
inoise = cv2.GaussianBlur(gray,(3,3),sigmaX=0)

#Canny
lowThresh = 50
upThresh = 200
edge = cv2.Canny(inoise,lowThresh,upThresh)

#Hough Line Transform
lines = cv2.HoughLines(edge,rho=1,theta=np.pi/180,threshold=120)
#rho = resolution of r, here 1 pixel
#theta = resolution of theta, here 1 degree
#threshold =  The minimum number of intersections to detect a line

img2 = img.copy()
#Plot detected lines
for r,theta in lines[0]:
    #(x0,y0) at the intersection point
    x0 = r*np.cos(theta)
    y0 = r*np.sin(theta)
    #(x1,y1) along the line to the left of (x0,y0) for 1000 units
    x1 =  int(x0-1000*np.sin(theta))
    y1 =  int(y0+1000*np.cos(theta))
    #(x2,y2) along the line to the right of (x0,y0) for 1000 units
    x2 =  int(x0+1000*np.sin(theta))
    y2 =  int(y0-1000*np.cos(theta))
    cv2.line(img2,(x1,y1),(x2,y2),(0,255,0),1)

cv2.imshow("Original",img)
cv2.imshow("Edge",edge)
cv2.imshow("Hough Lines",img2)
cv2.waitKey()
cv2.destroyAllWindows()

ให้สังเกตลูปในการพลอตเส้นตรง

ปกติเราสามารถกำหนดจุดสองจุดในการพลอตเส้นตรงใดๆ ได้

เช่นจากรูปข้างต้น เราอาจกำหนดจุดในการพลอตคือ
x1 = 0, y1 = r/sin(T)
x2 = r/cos(T), y2 =0
เพื่อเป็นการพลอตเส้นตรงจากขอบรูปไปยังอีกขอบหนึ่ง
แต่ปัญหาที่อาจจะเกิดขึ้นก็คือ ถ้าเป็นเส้นตรงในแนวนอน (T=0) ค่า y1 จะคำนวณไม่ได้

ดังนั้นในโค้ดของโปรแกรม จึงใช้อีกรูปแบบหนึ่ง ดังรูป

นั่นคือ
x1 = x0 - 1000 sin(T)
y1 = y0 + 1000 cos(T)

x2 = x0 + 1000 sin(T)
y2 = y0 - 1000 cos(T)

โดยที่ค่า 1000 เป็นระยะโดยประมาณครับ สามารถปรับเปลี่ยนได้

ข้อสังเกต

1. ความแม่นยำขึ้นอยู่กับวิธีการหาขอบ (edge detection)
2. ขอบเส้นตรงที่หาเจอ อาจจะไม่ใช่ขอบที่ถูกต้อง เพราะวิธีการนี้ตรวจสอบแค่ว่ามีพิกเซลที่วางตัวอยู่ในแนวเส้นตรงเดียวกันมากน้อยแค่ไหน ดังนั้น ในแนวใดๆ แม้ว่าจะไม่ใช่ขอบเส้นตรงจริงๆ หากมีจำนวนพิกเซลที่เป็นขอบอยู่จำนวนมาก ก็จะถูกกำหนดให้เป็นขอบเส้นตรง ตัวอย่างเช่น ถ้าเป็นต้นไม้พุ่มที่มีใบไม้จำนวนมาก ขอบของใบไม้หลายๆใบที่อยู่ในแนวเดียวกัน อาจจะทำให้ถูกคิดเป็นขอบเส้นตรงได้ ทั้งๆที่ไม่ใช่ เช่น เส้นสีเขียวในรูปด้านล่าง
3. การแก้ปัญหาในข้อ 2 อาจทำได้ เช่น กำหนดระยะห่างสูงสุดระหว่างเส้นสองเส้นในแนวเดียวกัน ถ้าระยะห่างนี้น้อยกว่าที่กำหนดก็ให้ถือว่าเส้นทั้งสองเป็นเส้นตรงเส้นเดียวกัน

การหาขอบ (Edge detection)

การหาขอบของภาพ ถูกนำมาใช้เพื่อหลายวัตถุประสงค์ เช่น เพื่อ segment ภาพนั้นออกเป็นส่วนๆ หรือ เพื่อเป็นกระบวนการเบื้องต้นสำหรับการหาเส้นตรง หรือวงกลมในภาพ ฯลฯ ในที่นี้เราจะลองใช้คำสั่งสำเร็จรูปของ OpenCV เพื่อหาขอบภาพ ด้วยสามวิธี ได้แก่

Canny

import cv2
import numpy as np

img = cv2.imread("lena.jpg",cv2.CV_LOAD_IMAGE_GRAYSCALE)
#Appy Gaussian blur to remove some noises
inoise = cv2.GaussianBlur(img,(3,3),sigmaX=0)

#Canny edge detection
lowThresh = 50
upThresh = 100
out = cv2.Canny(inoise,lowThresh,upThresh)

cv2.imshow("Origin",img)
cv2.imshow("Edge Detection",out)
cv2.waitKey()
cv2.destroyAllWindows()

Sobel

#Gradient X
gradX = cv2.Sobel(inoise,cv2.CV_16S,dx=1,dy=0)
#Gradient Y
gradY = cv2.Sobel(inoise,cv2.CV_16S,dx=0,dy=1)
#Convert Gradients to 8 bits
gradX = cv2.convertScaleAbs(gradX)
gradY = cv2.convertScaleAbs(gradY)
#Total Gradient (approximate)
#grad = abs(gradX) + abs(gradY)
out = cv2.addWeighted(gradX,1,gradY,1,0)

Laplacian

out = cv2.Laplacian(inoise,ddepth=cv2.CV_16S,ksize=3)
#Need to set depth to 16-bit signed integer because Laplacian can give negative intensity.
# Convert it to 8-bit image with absolute value.
out = cv2.convertScaleAbs(out)

Filtering / Convolution

การทำ Filtering / Convolution ก็คือการเอา kernel / template / filter (หน้าต่างขนาดเล็กว่ารูป มักมีขนาดเป็นเลขคี่ เช่น 3x3, 5x5, 7x7 ฯลฯ) มาคูณกับพื้นที่ของรูปที่ kernel นั้นซ้อนทับอยู่ แล้วหาผลรวม จากนั้นก็นำไปแทนที่พิกเซลที่อยู่ตรงกับตำแหน่งกลางของ kernel นั้น

เราลองมาดูการทำ filtering ด้วย average filter เพื่อเบลอรูปกันครับ

import cv2
import numpy as np

img = cv2.imread("lena.jpg")

ksize = 5   #kernel size
#create average kernel
kernel = np.ones((ksize,ksize))/(ksize*ksize)
ddepth = -1     #depth of output image, -1 is same as input image
out = cv2.filter2D(img,-1,kernel)
#same as: out = cv2.blur(img,(ksize,ksize))

cv2.imshow("Origin",img)
cv2.imshow("Result",out)
cv2.waitKey()
cv2.destroyAllWindows()

ถ้าต้องการเขียนโค้ดเอง สมมติว่าเป็นรูปสีเทา และคำนวณเฉพาะพิกเซลที่ตำแหน่งของ kernel ยังอยู่ในรูป

import cv2
import numpy as np

img = cv2.imread("lena.jpg", cv2.CV_LOAD_IMAGE_GRAYSCALE)

ksize = 5   #kernel size, odd value
#create average kernel
kernel = np.ones((ksize,ksize))/(ksize*ksize)

hksize = ksize%2    #half of kernel size
row,column = img.shape[:2]
out = img.copy()
for r in range(row-ksize+1):
    for c in range(column-ksize+1):
        temp = kernel*img[r:r+ksize,c:c+ksize]
        pix = temp.sum()
        out[r+hksize,c+hksize] = np.uint8(pix)

cv2.imshow("Origin",img)
cv2.imshow("Result",out)
cv2.waitKey()
cv2.destroyAllWindows()

ต้องระวังบรรทัด out[r+hksize,c+hksize] = np.uint8(pix) ให้ดีครับ เพราะถ้าไม่ใช่ average filter ในบรรทัดนี้อาจจะต้องทำการ clip/normalize ค่าพิกเซลให้อยู่ในช่วง 0-255 ด้วย

ถ้าเป็นรูปสี BGR ก็ต้องเพิ่มอีก 1 loop สำหรับ channel

hksize = ksize%2    #half of kernel size
row,column,channel = img.shape
out = img.copy()
for k in range(channel):
    for r in range(row-ksize+1):
        for c in range(column-ksize+1):
            temp = kernel*img[r:r+ksize,c:c+ksize,k]
            pix = temp.sum()
            out[r+hksize,c+hksize,k] = np.uint8(pix)

วิธีนี้จะใช้เวลาค่อนข้างมากครับ แต่ทำให้เราสามารถประยุกต์ใช้ kernel ตามที่เราต้องการได้

Thresholding การเปลี่ยนรูปเทาให้เป็นขาวดำ

ทำได้หลายวิธีครับ ในที่นี้เราจะใช้วิธีง่ายๆโดยการกำหนดค่า threshold ซึ่งถ้าค่าพิกเซลมากกว่าหรือเท่ากับค่านี้ก็จะเปลี่ยนให้เป็นสีขาว (255) มิฉะนั้นก็เปลี่ยนเป็นสีดำ (0)

import cv2

img = cv2.imread("lena.jpg",cv2.CV_LOAD_IMAGE_GRAYSCALE)
thres = 128
ret,bw = cv2.threshold(img,thres,255,cv2.THRESH_BINARY)

cv2.imshow("Original",img)
cv2.imshow("Binary",bw)

cv2.waitKey()
cv2.destroyAllWindows()

หรือจะใช้เทคนิคของ NumPy โดยเปลี่ยนโค้ดข้างต้นเป็น

import cv2

img = cv2.imread("lena.jpg",cv2.CV_LOAD_IMAGE_GRAYSCALE)
thres = 128
bw = 255*np.ones_like(img)
bw[img < thres] = 0

cv2.imshow("Original",img)
cv2.imshow("Binary",bw)

cv2.waitKey()
cv2.destroyAllWindows()

จะได้ผลลัพธ์เหมือนกันครับ