widzenie-komputerowe-MP/wko-05.ipynb

Widzenie komputerowe

05. Transformacje geometryczne i cechy obrazów [laboratoria]

Andrzej Wójtowicz (2021)

W poniższych materiałach zobaczymy jak przy pomocy biblitoeki OpenCV realizować transformacje geometryczne obrazu, wyszukiwać "ciekawe" elementy obrazu oraz jak łączyć je z innymi podobnymi elementami na innych obrazach.

Na początku załadujmy niezbędne biblioteki.

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Transformacje obrazu

Zacznijmy od podstawowych przekształceń geometrycznych. Aby efekt poszczególnych operacji był widoczny, dodamy pustą przestrzeń wokół testowego obrazu, na którym będziemy pracować:

lena = cv.imread("img/lena.png", cv.IMREAD_COLOR)
image = np.zeros((1024, 1024, 3), np.uint8)

for i in range(3):
    image[256:768, 256:768, i] = lena[:,:, i];

plt.figure(figsize=(5,5))
plt.imshow(image[:,:,::-1]);

W OpenCV transformacje afiniczne są wykonywane przy pomocy definicji macierzy 2x3 (zamiast bardziej klasycznej 3x3), która przekazywana jest do funkcji cv.warpAffine(). Poniżej mamy przykład translacji (przesunięcia) o 100 pikseli w prawo wzdłuż osi _x i o 150 pikseli wzdłuż osi y:

mat = np.float32([
    [1.0, 0.0, 100],
    [0.0, 1.0, 150]
])

image_translated = cv.warpAffine(image, mat, image.shape[0:2])

plt.figure(figsize=(5,5))
plt.imshow(image_translated[:,:,::-1]);

Zobaczmy co się stanie po przeskalowaniu dwukrotnym wzdłuż osi _x:

mat = np.float32([
    [2.0, 0.0, 0],
    [0.0, 1.0, 0]
])

image_scaled = cv.warpAffine(image, mat, image.shape[0:2])

plt.figure(figsize=(5,5))
plt.imshow(image_scaled[:,:,::-1]);

Nasz wyjściowy obraz "wyszedł" poza wymiary wejściowe, więc musimy podać w funkcji trochę większe rozmiary wynikowego obrazu (w tym wypadku chodzi o szerokość):

mat = np.float32([
    [2.0, 0.0, 0],
    [0.0, 1.0, 0]
])

image_scaled = cv.warpAffine(image, mat, (2*image.shape[1], image.shape[0]))

plt.figure(figsize=(10,10))
plt.imshow(image_scaled[:,:,::-1]);

Poniżej mamy dwukrotne przeskalowanie w obu kierunkach (efekt jest zasadniczo widoczny patrząc na wynikowe skale obu osi):

mat = np.float32([
    [2.0, 0.0, 0],
    [0.0, 2.0, 0]
])

image_scaled = cv.warpAffine(image, mat, (2*image.shape[1], 2*image.shape[0]))

plt.figure(figsize=(5,5))
plt.imshow(image_scaled[:,:,::-1]);

W przypadku obrotu musimy pamiętać, że mamy inaczej zdefiniowany układ współrzędnych (oś _y), więc uzyskujemy trochę inną formę macierzy obrotu:

angle = 30 * np.pi / 180.0

cos_theta = np.cos(angle)
sin_theta = np.sin(angle)

mat = np.float32([
    [ cos_theta, sin_theta, 0],
    [-sin_theta, cos_theta, 0]
])

image_rotated = cv.warpAffine(image, mat, image.shape[0:2])

plt.figure(figsize=(5,5))
plt.imshow(image_rotated[:,:,::-1]);

Powyżej mieliśmy obrót wokół środka układu, ale możemy też uzyskać obrót wokół wskazanego punktu, np. środka obrazu:

angle = 30 * np.pi / 180.0

cos_theta = np.cos(angle)
sin_theta = np.sin(angle)

center_x = image.shape[0] / 2
center_y = image.shape[1] / 2

t_x = (1 - cos_theta) * center_x - sin_theta * center_y
t_y = sin_theta * center_x + (1 - cos_theta) * center_y

mat = np.float32([
    [ cos_theta, sin_theta, t_x],
    [-sin_theta, cos_theta, t_y]
])

image_rotated = cv.warpAffine(image, mat, image.shape[0:2])

plt.figure(figsize=(5,5))
plt.imshow(image_rotated[:,:,::-1]);

Na szczęście zamiast samodzielnego liczenia macierzy obrotu można użyć funkcji cv.getRotationMatrix2D() (możemy też jej użyć do skalowania):

mat = cv.getRotationMatrix2D((center_x, center_y), 30, 1)

image_rotated = cv.warpAffine(image, mat, image.shape[0:2])

plt.figure(figsize=(5,5))
plt.imshow(image_rotated[:,:,::-1]);

Poniżej mamy przykład pochylenia:

mat = np.float32([
    [1.0, 0.1, 0],
    [0.0, 1.0, 0]
])

image_sheared = cv.warpAffine(image, mat, image.shape[0:2])

plt.figure(figsize=(5,5))
plt.imshow(image_sheared[:,:,::-1]);

Jeżeli mamy kilka transformacji do wykonania, to bardziej efektywnym rozwiązaniem będzie uprzednie wymnożenie macierzy odpowiedzialnych za poszczególne transformacje (zamiast sukcesywnego wykonywania pojedynczych, izolowanych operacji):

m_scale = np.float32([
    [0.5, 0.0],
    [0.0, 1.0]
])

m_shear = np.float32([
    [1.0, 0.1],
    [0.0, 1.0]
])

m_rotation = np.float32([
    [ cos_theta, sin_theta],
    [-sin_theta, cos_theta]
])

v_translation = np.float32([
    [150],
    [90]
])

mat_tmp = m_rotation @ m_shear @ m_scale

mat_transformation = np.append(mat_tmp, v_translation, 1)

image_transformed = cv.warpAffine(image, mat_transformation, image.shape[0:2])

plt.figure(figsize=(5,5))
plt.imshow(image_transformed[:,:,::-1]);

Dla punktów o konkretnych współrzędnych możemy też obliczyć ich docelowe współrzędne po transformacjach:

src_points = np.float32([[50, 50], [50, 249], [249, 50], [249, 249]])
dst_points = (mat_tmp @ src_points.T + v_translation).T

print(dst_points)

[[200.98076 118.30127]
 [317.71466 280.6903 ]
 [287.15027  68.55127]
 [403.88422 230.94032]]

Mając co najmniej 3 punkty wejściowe i ich odpowiedniki wyjściowe, możemy oszacować macierz transformacji przy pomocy funkcji cv.estimateAffine2D():

mat_estimated_1 = cv.estimateAffine2D(src_points[:3], dst_points[:3])[0]
print("Explicit matrix:\n\n", mat_transformation)
print("\nEstimated matrix:\n\n", mat_estimated_1)

Explicit matrix:

 [[  0.4330127    0.58660257 150.        ]
 [ -0.25         0.8160254   90.        ]]

Estimated matrix:

 [[  0.43301261   0.58660252 150.00000192]
 [ -0.25         0.81602532  90.00000368]]

Mając więcej punktów możemy spróbować uzyskać dokładniejsze oszacowanie:

mat_estimated_2 = cv.estimateAffine2D(src_points, dst_points)[0]
print("Explicit matrix:\n\n", mat_transformation)
print("\nEstimated matrix:\n\n", mat_estimated_2)

Explicit matrix:

 [[  0.4330127    0.58660257 150.        ]
 [ -0.25         0.8160254   90.        ]]

Estimated matrix:

 [[  0.43301273   0.58660264 149.99997897]
 [ -0.24999996   0.81602536  89.99999603]]

Załóżmy, że chcielibyśmy teraz uzyskać zmianę perspektywy, co w efekcie spowodowałoby, że nasz kwadratowy obraz zamieniłby się w trapez:

image_to_transform = np.zeros(image.shape, dtype=np.float32)
dst_points = np.float32([[356, 256], [667, 256], [767, 767], [256, 767]])
cv.fillConvexPoly(image_to_transform, np.int32(dst_points), (0.8, 0.8, 0.5), cv.LINE_AA);

plt.figure(figsize=[10,10])
plt.subplot(121)
plt.imshow(image[:,:,::-1])
plt.title('Original image')

plt.subplot(122)
plt.imshow(image_to_transform[:,:,::-1])
plt.title('Destination plane');

Niestety przy pomocy transformacji afinicznych nie jesteśmy w stanie uzyskać zakładanego efektu:

src_points = np.float32([[256, 256], [767, 256], [767, 767], [256, 767]])
mat_estimated = cv.estimateAffine2D(src_points, dst_points)[0]

image_transformed_a = cv.warpAffine(image, mat_estimated, image.shape[0:2])

plt.figure(figsize=[15,10])
plt.subplot(131)
plt.imshow(image[:,:,::-1])
plt.title('Original image')
plt.subplot(132)
plt.imshow(image_to_transform[:,:,::-1])
plt.title('Destination plane')
plt.subplot(133)
plt.imshow(image_transformed_a[:,:,::-1])
plt.title('Transformed image (affine)');

Do uzyskania zakładanego efektu musimy znaleźć macierz homografii 3x3 przy pomocy co najmniej 4 punktów i funkcji cv.findHomography():

mat_h = cv.findHomography(src_points, dst_points)[0]
print(mat_h)

[[ 5.08838663e-01 -3.27547619e-01  2.51229024e+02]
 [ 0.00000000e+00  3.44904761e-01  1.25737302e+02]
 [ 0.00000000e+00 -6.40366802e-04  1.00000000e+00]]

Przy pomocy funkcji cv.warpPerspective() możemy dokonać teraz zmiany perspektywy:

image_transformed_h = cv.warpPerspective(image, mat_h, image.shape[0:2])

plt.figure(figsize=[20,10])
plt.subplot(141)
plt.imshow(image[:,:,::-1])
plt.title('Original image')
plt.subplot(142)
plt.imshow(image_to_transform[:,:,::-1])
plt.title('Destination plane')
plt.subplot(143)
plt.imshow(image_transformed_a[:,:,::-1])
plt.title('Transformed image (affine)')
plt.subplot(144)
plt.imshow(image_transformed_h[:,:,::-1])
plt.title('Transformed image (homography)');

Poprzez znalezienie macierzy homografii i zmianę perspektywy możemy np. wyrównywać zdjęcia dokumentów w celu ich dalszej analizy. Poniżej mamy przykład z oznaczeniem okładki zeszytu:

book = cv.imread("img/caribou.jpg", cv.IMREAD_COLOR)
image_h, image_w = book.shape[0:2]

src_points = np.array([[210, 190], [380, 280], [290, 500], [85, 390]], dtype=float)
dst_points = np.array([[0, 0],[image_w-1, 0], [image_w-1, image_h-1], [0, image_h-1]], dtype=float)

colors = [(0,0,255), (0,255,0), (255,0,0), (90,200,200)]
for (x,y), i in(zip(src_points, colors)):
    book = cv.circle(book, (int(x),int(y)), 10, i, -1)

mat = cv.findHomography(src_points, dst_points)[0]

book_aligned = cv.warpPerspective(book, mat, (image_w, image_h))

plt.figure(figsize=[10,10])
plt.subplot(121)
plt.imshow(book[:,:,::-1])
plt.title('Original image')
plt.subplot(122)
plt.imshow(book_aligned[:,:,::-1])
plt.title('Aligned image');

W przypadku dokumentów, możemy np. dysponować wzorem dokumentu do uzupełnienia oraz zdjęciem wypełnionego dokumentu. Po wyrówaniu zdjęcia moglibyśmy łatwiej analizować wypełnione pola:

Cztery niezbędne punkty moglibyśmy wybierać ręcznie, jednak w przypadku masowego przetwarzania lepszym rozwiązaniem jest automatyczne znalezienie odpowiadających sobie punktów (niekoniecznie rogów dokumentu) na wzorcu i przetwarzanym zdjęciu. Poniżej zobaczymy jak można znaleźć i dopasować takie punkty.

Punkty kluczowe

Punkty kluczowe (ang. _keypoints) są punktami na obrazie, o których możemy myśleć w kategoriach dobrze zdefiniowanych narożników, relatywnie odpornych na przekształcenia obrazu. Każdy punkt kluczowy ma przypisany descryptor tj. wektor cech.

Wczytajmy przykładowy obraz, na którym postaramy się znaleźć takie punkty kluczowe:

python_book = cv.imread("img/book-python-cover.jpg", cv.IMREAD_COLOR)
python_book_gray = cv.cvtColor(python_book, cv.COLOR_BGR2GRAY)

plt.figure(figsize=[10,10])
plt.subplot(121)
plt.imshow(python_book[:,:,::-1])
plt.title('Original image')
plt.subplot(122)
plt.imshow(python_book_gray, cmap='gray')
plt.title('Grayscale image');

Jednym z najpopularniejszych (i nieopatentowanych) algorytmów wykrywania i opisywania punktów kluczowych jest ORB (dokumentacja). Poniżej znajduje się wizualizacja wykrytych punktów kluczowych:

orb_alg = cv.ORB_create()

keypoints_1 = orb_alg.detect(python_book_gray, None)

keypoints_1, descriptors_1 = orb_alg.compute(python_book_gray, keypoints_1)

python_book_keypoints = cv.drawKeypoints(python_book, keypoints_1, None, color=(0,255,0), 
                                         flags=cv.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)

plt.figure(figsize=(10,10))
plt.imshow(python_book_keypoints[:,:,::-1])
plt.title('Image keypoints');

Łączenie punktów kluczowych

Załóżmy, że chcemy sprawdzić czy na danym zdjęciu znajduje się okładka powyższej książki:

book_in_hands = cv.imread("img/book-python-in-hands.jpg", cv.IMREAD_COLOR)
book_in_hands_gray = cv.cvtColor(book_in_hands, cv.COLOR_BGR2GRAY)

plt.figure(figsize=[10,10])
plt.subplot(121)
plt.imshow(book_in_hands[:,:,::-1])
plt.title('Original image')
plt.subplot(122)
plt.imshow(book_in_hands_gray, cmap='gray')
plt.title('Grayscale image');

Na początku znajdźmy i opiszmy punkty kluczowe. Oba kroki możemy wykonać jednocześnie przy pomocy metody detectAndCompute():

orb_alg = cv.ORB_create(nfeatures=1000)

keypoints_1, descriptors_1 = orb_alg.detectAndCompute(python_book_gray, None)
keypoints_2, descriptors_2 = orb_alg.detectAndCompute(book_in_hands_gray, None)

Jeśli uda nam się połączyć odpowiadające sobie punkty kluczowe i wyliczyć macierz homografii, to z dość dużą dozą pewności będziemy mogli przyjąć, że na zdjęciu znajduje się szukany obiekt i będziemy mogli z grubsza wskazać gdzie on się znajduje. Punkty kluczowe są opisane przez wektory, więc np. moglibyśmy _brute-forcem i odległością Hamminga poszukać najlepszych odpowiedników. W poniższym przypadku pozostawimy 5% najlepszych dopasowań i zwizualizujemy te dopasowania:

matcher = cv.DescriptorMatcher_create(cv.DESCRIPTOR_MATCHER_BRUTEFORCE_HAMMING)
matches = matcher.match(descriptors_1, descriptors_2, None)
matches = list(matches)
matches.sort(key=lambda x: x.distance, reverse=False) # sort by distance

num_good_matches = int(len(matches) * 0.05) # save 5%
matches = matches[:num_good_matches]

image_matched = cv.drawMatches(python_book, keypoints_1, book_in_hands, keypoints_2, matches, None)
plt.figure(figsize=(10,10))
plt.imshow(image_matched[...,::-1]);

Jak widać część punktów jest niepoprawnie dopasowana. Na szczęście obliczenie macierzy homografii odbywa się przy pomocy algorytmu RANSAC (ang. _random sample consensus), co pozwala wziąć pod uwagę, że pewne odzwierciedlenia są złe. Aby oznaczyć miejsce występowania obiektu, transformujemy współrzędne narożników źródłowego zdjęcia przy użyciu otrzymanej macierzy i funkcji cv.perspectiveTransform():

src_pts = np.float32([keypoints_1[m.queryIdx].pt for m in matches]).reshape(-1,1,2)
dst_pts = np.float32([keypoints_2[m.trainIdx].pt for m in matches]).reshape(-1,1,2)

mat, mask = cv.findHomography(src_pts, dst_pts, cv.RANSAC, 5.0)
matches_mask = mask.ravel().tolist()
height, width, _ = python_book.shape

pts = np.float32([[0,0], [0,height-1], [width-1,height-1], [width-1,0]]).reshape(-1,1,2)

dst = cv.perspectiveTransform(pts, mat)

book_in_hands_poly = cv.polylines(book_in_hands,
                                  [np.int32(dst)], True, (0,0,255), 
                                  10, cv.LINE_AA)
    
draw_params = dict(matchColor=(0,255,0),
                   matchesMask=matches_mask, # draw only inliers
                   flags=2) # NOT_DRAW_SINGLE_POINTS
image_matched = cv.drawMatches(python_book, keypoints_1, 
                               book_in_hands_poly, keypoints_2, 
                               matches, None, **draw_params)

plt.figure(figsize=(10,10))
plt.imshow(image_matched[...,::-1]);

Zadanie 1

Napisz interaktywną aplikację przy pomocy HighGUI, która pozwoli na obrazie img/billboards.jpg wybrać 4 punkty i umieścić w zaznaczonym obszarze wybrany obraz, np. img/bakery.jpg.

Zadanie 2

Załóżmy, że mamy dwa zdjęcia, które chcielibyśmy połączyć w panoramę:

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
boat_1 = cv.imread("img/boat_1.jpg", cv.IMREAD_COLOR)
boat_2 = cv.imread("img/boat_2.jpg", cv.IMREAD_COLOR)

plt.figure(figsize=[10,10])
plt.subplot(121)
plt.imshow(boat_1[:,:,::-1])
plt.title('Image 1')
plt.subplot(122)
plt.imshow(boat_2[:,:,::-1])
plt.title('Image 2');

Standardowo w OpenCV panoramę tworzy się przy pomocy klasy Stitcher, której można podać listę zdjęć do połączenia:

stitcher_alg = cv.Stitcher_create()
_, panorama = stitcher_alg.stitch([boat_1, boat_2])
plt.figure(figsize=[20,10]) 
plt.imshow(panorama[:,:,::-1])
plt.title('Panorama');

Spróbuj zaimplementować samodzielne stworzenie panoramy. Poniższy efekt będzie wystarczający, aczkolwiek możesz również spróbować zniwelować widoczną różnicę w intensywnościach na granicy zdjęć.

class Image_Stitching():
    def __init__(self) :
        self.ratio=0.85
        self.min_match=10
        self.sift=cv.xfeatures2d.SIFT_create()
        self.smoothing_window_size=800

    def registration(self,img1,img2):
        kp1, des1 = self.sift.detectAndCompute(img1, None)
        kp2, des2 = self.sift.detectAndCompute(img2, None)
        matcher = cv.BFMatcher()
        raw_matches = matcher.knnMatch(des1, des2, k=2)
        good_points = []
        good_matches=[]
        for m1, m2 in raw_matches:
            if m1.distance < self.ratio * m2.distance:
                good_points.append((m1.trainIdx, m1.queryIdx))
                good_matches.append([m1])
        img3 = cv.drawMatchesKnn(img1, kp1, img2, kp2, good_matches, None, flags=2)
        cv.imwrite('matching.jpg', img3)
        if len(good_points) > self.min_match:
            image1_kp = np.float32(
                [kp1[i].pt for (_, i) in good_points])
            image2_kp = np.float32(
                [kp2[i].pt for (i, _) in good_points])
            H, status = cv.findHomography(image2_kp, image1_kp, cv.RANSAC,5.0)
        return H

    def create_mask(self,img1,img2,version):
        height_img1 = img1.shape[0]
        width_img1 = img1.shape[1]
        width_img2 = img2.shape[1]
        height_panorama = height_img1
        width_panorama = width_img1 +width_img2
        offset = int(self.smoothing_window_size / 2)
        barrier = img1.shape[1] - int(self.smoothing_window_size / 2)
        mask = np.zeros((height_panorama, width_panorama))
        if version== 'left_image':
            mask[:, barrier - offset:barrier + offset ] = np.tile(np.linspace(1, 0, 2 * offset ).T, (height_panorama, 1))
            mask[:, :barrier - offset] = 1
        else:
            mask[:, barrier - offset :barrier + offset ] = np.tile(np.linspace(0, 1, 2 * offset ).T, (height_panorama, 1))
            mask[:, barrier + offset:] = 1
        return cv.merge([mask, mask, mask])

    def blending(self,img1,img2):
        H = self.registration(img1,img2)
        height_img1 = img1.shape[0]
        width_img1 = img1.shape[1]
        width_img2 = img2.shape[1]
        height_panorama = height_img1
        width_panorama = width_img1 +width_img2

        panorama1 = np.zeros((height_panorama, width_panorama, 3))
        mask1 = self.create_mask(img1,img2,version='left_image')
        panorama1[0:img1.shape[0], 0:img1.shape[1], :] = img1
        panorama1 *= mask1
        mask2 = self.create_mask(img1,img2,version='right_image')
        panorama2 = cv.warpPerspective(img2, H, (width_panorama, height_panorama))*mask2
        result=panorama1+panorama2

        rows, cols = np.where(result[:, :, 0] != 0)
        min_row, max_row = min(rows), max(rows) + 1
        min_col, max_col = min(cols), max(cols) + 1
        final_result = result[min_row:max_row, min_col:max_col, :]
        return final_result
    
def main(img1,img2):
    final=Image_Stitching().blending(img1,img2)
    cv.imwrite('panorama.jpg', final)

main(boat_1, boat_2)
panorama_jpg = cv.imread("panorama.jpg", cv.IMREAD_COLOR)
plt.figure(figsize=[10,10])
plt.imshow(panorama_jpg[:,:,::-1])
plt.title('Panorama')
    

[ WARN:0@11.790] global /io/opencv_contrib/modules/xfeatures2d/misc/python/shadow_sift.hpp (13) SIFT_create DEPRECATED: cv.xfeatures2d.SIFT_create() is deprecated due SIFT tranfer to the main repository. https://github.com/opencv/opencv/issues/16736

Text(0.5, 1.0, 'Panorama')