homography matrix 다루기 / cv2.findHomography( ), cv2.warpPerspective( )

 

 

 제목을 뭐로 해야할지 몰라서 homography matrix 다루기라 두루뭉실하게 적어놨다. 오늘은 아래의 그림과 목차 흐름으로 정리하고자 한다.

 

1. cv2.xfeatures2d.SIFT_create( ) 함수로 keypoint와 descriptor 추출
2. cv2.BFMatcher( ) 함수로 keypoint 매칭
3. cv2.findHomography( ) 함수로 homography matrix, H 추출
4. 추출한 H로 이미지 변환
   → cv2.warpPerspective( ) 함수 사용하면 쉽게 할 수 있는데,
      이번 글에서는 H를 직접 다루는 방법으로 정리하고자 한다.
5. 정합 결과 확인

 

 먼저 말해두자면, 난 SIFT 이론에 대해 아직은 정확히 모른다. -> 20.01.15 논문 쓰려고 공부했다. 이해한 데까지 정리도 했다. [openCV] - SIFT 알고리즘 여기를 보면 된다!

 

 homography matrix 이론에 대해서도 정확하게 아는 건 아니지만 개발에 필요해서 구현 방법을 중심으로 짧게 나마 알아봤고, 그 내용들에 대해 정리해놓고자 한다.

 

 SIFT 사용법과 matching 방법은 앞선 글 cv2.BFMatcher( )에 대해 알아보자 에 설명해뒀으니 넘어가도록 하겠다.

 

 

---

 

 

1. Homography Matrix, H

 3x3 행렬로 변환 행렬에 해당되는 H는 아래와 같이 표현되며 cv2.findHomography( ) 함수를 통해 구해줄 수 있다.

 

 

 

 image A에서 뽑은 keypoint와 매칭되는 image B의 keypoint를 cv2.findHomography( ) 함수에 넣어주면 된다. 이렇게 구한 H를 사용하면 우리는 image A와 image B 정합시켜줄 수 있다. 더 정확히 표현하면 image A를 image B에 겹치게 변환시켜줄 수 있다. 설명이 헷갈릴텐데 밑에 코드와 결과를 같이 보면 충분히 이해할 수 있을 거다.

 

 변환 과정에는 H와 cv2.warpPerspective( ) 함수를 사용하면 되는데 이번 글에선 OpenCV 함수를 안 쓰고 직접 구현하는 방법으로 정리하려고 한다. ( cv2.warpPerspective( ) 함수 사용법은 구글링 하면 많이 나온다. )

 

 

 방법은 생각보다 간단하다. H와 좌표계를 행렬곱 연산을 해주면 된다. 위의 수식은 homography matrix를 검색해보면 자주 검색되는 식인데, 중요한 점은 위의 식 그대로 연산을 진행해주면 안 되고 아래 식으로 계산해야한다는 점이다.

 

 

x'를 정확하게 구하려면

h11*x    +    h12*y    +    h13 에서 계산을 끝낼 게 아니라

이를 h31*x    + h32*y    +    1로 나눠줘야 한다.

그래야 올바른 결과를 얻을 수 있다.

 

이론 설명은 이쯤에서 멈추고 코드를 통해 이해해보도록 하자.

 

 

 

---

 

 

< 왼쪽 사진 >

< 오른쪽 사진 >

import numpy as np
import cv2


imageA = cv2.imread('./pano_2.jpeg') # 오른쪽 사진
imageB = cv2.imread('./pano_1.jpeg') # 왼쪽 사진

grayA = cv2.cvtColor(imageA,cv2.COLOR_BGR2GRAY)
grayB = cv2.cvtColor(imageB,cv2.COLOR_BGR2GRAY)


sift = cv2.xfeatures2d.SIFT_create()
kpA, desA = sift.detectAndCompute(grayA, None)
kpB, desB = sift.detectAndCompute(grayB, None)


bf = cv2.BFMatcher()
matches = bf.match(desA, desB)


sorted_matches = sorted(matches, key = lambda x : x.distance)
res = cv2.drawMatches(imageA, kpA, imageB, kpB, sorted_matches[:30], None, flags = 2)


src = np.float32([kpA[m.queryIdx].pt for m in matches]).reshape((-1, 1, 2))
dst = np.float32([kpB[m.trainIdx].pt for m in matches]).reshape((-1, 1, 2))
H, status = cv2.findHomography(src, dst, cv2.RANSAC, 5.0)


before = []
for x in range(imageA.shape[1]):
    for y in range(imageA.shape[0]):
        point = [x, y, 1]
        before.append(point)
before = np.array(before).transpose()

after = np.matmul(H, before)
after = after / after[2, :]
after = after[:2, :]
after = np.round(after, 0).astype(np.int)


height, width, _ = imageA.shape
result = np.zeros((height, width * 2, 3), dtype = np.uint8)
for pt1, pt2 in zip(before[:2, :].transpose(), after.transpose()):
    if pt2[1] >= height:
        continue

    if np.sum(pt2 < 0) >= 1:
        continue
    
    result[pt2[1], pt2[0]] = imageA[pt1[1], pt1[0]]
result[0: height, 0 : width] = imageB

cv2.imshow('result', result)

• [Line 5, 6]
  우리는 왼쪽사진은 그대로 두고 오른쪽 사진을 왼쪽 사진과 잘 맞게 변형을 할 것이다.
  그래서 코드의 전체 흐름은 오른쪽 사진을 image A로 두고 image A를 변환하는 흐름으로 진행된다.
• [Line 8 ~ 22]
  앞선 글에서 다뤘던 내용이므로 설명은 생략하겠다.
• [Line 25 ~ 27]
  cv2.findHomography( )를 통해 H를 구하는 과정이다.
  image A ( = 오른쪽 사진 )의 src ( = image A의 keypoint 좌표 )는
  image B ( = 왼쪽 사진 )의 dst ( = image B의 keypoint 좌표 )와 매칭되므로
  src가 dst로 변환될 수 있도록 해주는 변환 행렬 H를 cv2.findHomography( ) 함수가 구해준다.
  ※ cv2.RANSAC, status 에 대해서는 다른 글에서 정리하도록 하겠다.
• [Line 30 ~ 35]
  변환 전 좌표들에 해당된다.
  의 모양에 맞춰주기 위한 과정이다.
• [Line 37 ~ 40]
  H와 변환 전 좌표 before를 행렬곱 연산을 통해 변환 후 좌표 after를 얻는 과정이다.
  Line 38은 위에서 말한  을 나눠주는 과정에 해당하고,
  Line 39는 에서 와 만을 뽑는 과정이다.
• [Line 43 ~ 53]
  최종적으로 변환 결과를 그려주는 과정이다.
  변환 후 좌표가 가능 범위를 벗어나는 경우에 에러가 날 수 있기 때문에 continue를 넣어줬다.

 

결과 이미지는 다음과 같다.

 

 

 왼쪽 이미지에 오른쪽 이미지가 맞춰서 변형된 걸 확인할 수 있다. 검정색 빗살이 쳐지는 이유는 아직 모르겠다. cv2.warpPerspective( ) 함수를 이용해 정합하면 저런 현상이 안 나타나는 걸 보면 코드에서 살짝 부족한 점이 있는 것 같다. 보간법 사용 여부에 따라 결과 차이가 생기는 건가?

 

 

 

---

! 우선 광고 시간 !

728x90

 

---

 

 

2019.10.08. 내용 추가

 

 검정색 빗살이 쳐지는 현상 해결에 대해 정리한 블로그를 우연히 읽게 됐다. 위와 같은 현상이 발생하는 이유는, before point가 after point로 변환되면서 after image의 모든 pixel 영역을 채워주지 못하기 때문이라고 한다. 그리고 이와 같은 문제는 보간(interpolation)을 통해 해결할 수 있다. max interpolation 방법이 가장 쉽게 구현할 수 있을 것 같아 구현해봤다.

h, w, _ = result.shape
kernel_size = 5
inpterpolation = np.zeros((h, w, 3), dtype = np.uint8)
for x in range(int((w / kernel_size))):
    for y in range((int(h / kernel_size))):
        roi = result[int(y * kernel_size) : int(y * kernel_size) + kernel_size,
                     int(x * kernel_size) : int(x * kernel_size) + kernel_size]

        roi_max_b = roi[:,:,0].max()
        roi_max_g = roi[:,:,1].max()
        roi_max_r = roi[:,:,2].max()
        
        for i in range(kernel_size):
            for j in range(kernel_size):
                if np.sum(roi[i, j] == [0,0,0]) == 3:
                    roi[i, j] = [roi_max_b, roi_max_g, roi_max_r]
        
        inpterpolation[int(y * kernel_size) : int(y * kernel_size) + kernel_size,
                      int(x * kernel_size) : int(x * kernel_size) + kernel_size] = roi

 

• 각 채널 별로 max 값을 뽑아서 roi 이미지 중 [0, 0, 0]을 갖는 pixel에 대입하도록 구현했다.

 

 

 

 결과를 봐보면 검정색 빗살이 사라지긴 했는데 보간 결과가 그리 만족스럽진 않다. nearest 방식을 적용해보면 더 나으려나? 어떻게 구현하는 지 알게 되면 그때 또 추가로 정리해야겠다.

 

 

 

---

 

 

2020.07.10. 내용 추가

 

 마이클플리님이 댓글을 통해 남겨주신 조언을 토대로 알아본 결과 지금까지의 방식은 'Forward Mapping'이었고, 'Backward Mapping'을 통해 변환을 해주면 빗살 무늬 현상을 해결할 수 있다고 한다. 

 

[출처] https://www.researchgate.net/publication/323796200_Self-Supervised_Monocular_Image_Depth_Learning_and_Confidence_Estimation

 

 자세한 내용은 이 블로그를 참고하자. Backward Mapping은 다음과 같이 구현했다.

# backward mapping
before2 = []
for x in range(imageA.shape[1], imageA.shape[1] * 2):
    for y in range(imageA.shape[0]):
        point = [x, y, 1]
        before2.append(point)
before2 = np.array(before2).transpose()

Hinv = np.linalg.inv(H)

after2 = np.matmul(Hinv, before2)
after2 = after2 / after2[2, :]
after2 = after2[:2, :]
after2 = np.round(after2, 0).astype(np.int)

height, width, _ = imageA.shape
result2 = np.zeros((height, width * 2, 3), dtype = np.uint8)
for pt1, pt2 in zip(before2[:2, :].transpose(), after2.transpose()):
    if pt2[1] >= height or pt2[0] >= width:
        continue

    if np.sum(pt2 < 0) >= 1:
        continue
    
    result2[pt1[1], pt1[0]] = imageA[pt2[1], pt2[0]]
result2[0: height, 0 : width] = imageB
cv2.imshow('result2', result2)


# cv2.warpPerspective( )
result3 = cv2.warpPerspective(imageA, H, (imageA.shape[1] + imageB.shape[1], imageA.shape[0]))
result3[0 : imageA.shape[0], 0 : imageB.shape[1]] = imageB
cv2.imshow('result3', result3)

 

• [Line 9]
  변환된 이미지의 좌표계를 기준으로 변환 전 이미지의 픽셀값을 가져오는 흐름이라
  Homography 행렬의 역행렬을 구해 변환에 사용하였다.
• [Line 30 ~ 33]
  cv2.warpPespective( ) 함수를 통해 구한 결과 영상과 비교를 하기 위해 추가했다.

 

 

우선 Backward Mapping 결과 영상(result2)이다.

 

 

다음은 cv2.warpPespective( )를 사용했을 때의 결과 영상(result3)이다.

 

 

 Result2 영상이 보간을 사용했을 때 보다 훨씬 더 자연스러운 결과라는 것을 확인할 수 있다. 하지만 Result3에 비하면 살짝 부족한 느낌은 든다. Backward Mapping 단계에서 발생하는 소수형태의 좌표 결과를 어떻게 처리하느냐에 따라 영상의 퀄리티가 달라지는 것 같다. 나는 깊은 고민없이 단순히 반올림만 해준 것이라 아무래도 퀄리티가 상대적으로 낮은 것 같다. 나중에.. 좀 더 연구가 필요하다 생각이 들거나 우연히 또 알게된다면 시도해보독 하겠다.

 

 

 

[ 참고 사이트 ]

https://darkpgmr.tistory.com/79

https://b.mytears.org/2007/10/599/

https://wewinserv.tistory.com/89

https://yongchul-note.tistory.com/4

https://www.researchgate.net/publication/323796200_Self-Supervised_Monocular_Image_Depth_Learning_and_Confidence_Estimation