from PIL import Image

img = Image.open('lena.png')
img2 = img.rotate(45)       # 逆時(shí)針旋轉(zhuǎn)45°
img2.save("lena_rot45.png")
img2.show()

"""
圖像匹配——SIFT點(diǎn)特征匹配實(shí)現(xiàn)步驟：
    （1）讀取圖像；
    （2）定義sift算子；
    （3）通過sift算子對(duì)需要匹配的圖像進(jìn)行特征點(diǎn)獲?。?
        a.可獲取各匹配圖像經(jīng)過sift算子的特征點(diǎn)數(shù)目
    （4）可視化特征點(diǎn)（在原圖中標(biāo)記為圓圈）；
        a.為方便觀察，可將匹配圖像橫向拼接
    （5）圖像匹配（特征點(diǎn)匹配）；
        a.通過調(diào)整ratio獲取需要進(jìn)行圖像匹配的特征點(diǎn)數(shù)量（ratio值越大，匹配的線條越密集，但錯(cuò)誤匹配點(diǎn)也會(huì)增多）
        b.通過索引ratio選擇固定的特征點(diǎn)進(jìn)行圖像匹配
    （6）將待匹配圖像通過旋轉(zhuǎn)、變換等方式將其與目標(biāo)圖像對(duì)齊
"""

import cv2              # opencv版本需為3.4.2.16
import numpy as np      # 矩陣運(yùn)算庫
import time             # 時(shí)間庫

original_lena = cv2.imread('lena.png')          # 讀取lena原圖
lena_rot45 = cv2.imread('lena_rot45.png')       # 讀取lena旋轉(zhuǎn)45°圖

sift = cv2.xfeatures2d.SIFT_create()

# 獲取各個(gè)圖像的特征點(diǎn)及sift特征向量
# 返回值kp包含sift特征的方向、位置、大小等信息；des的shape為（sift_num， 128）， sift_num表示圖像檢測(cè)到的sift特征數(shù)量
(kp1, des1) = sift.detectAndCompute(original_lena, None)
(kp2, des2) = sift.detectAndCompute(lena_rot45, None)

# 特征點(diǎn)數(shù)目顯示
print("=========================================")
print("=========================================")
print('lena 原圖  特征點(diǎn)數(shù)目：', des1.shape[0])
print('lena 旋轉(zhuǎn)圖 特征點(diǎn)數(shù)目：', des2.shape[0])
print("=========================================")
print("=========================================")

# 舉例說明kp中的參數(shù)信息
for i in range(2):
    print("關(guān)鍵點(diǎn)", i)
    print("數(shù)據(jù)類型:", type(kp1[i]))
    print("關(guān)鍵點(diǎn)坐標(biāo):", kp1[i].pt)
    print("鄰域直徑:", kp1[i].size)
    print("方向:", kp1[i].angle)
    print("所在的圖像金字塔的組:", kp1[i].octave)

print("=========================================")
print("=========================================")
"""
首先對(duì)原圖和旋轉(zhuǎn)圖進(jìn)行特征匹配，即圖original_lena和圖lena_rot45
"""
# 繪制特征點(diǎn)，并顯示為紅色圓圈
sift_original_lena = cv2.drawKeypoints(original_lena, kp1, original_lena, color=(255, 0, 255))
sift_lena_rot45 = cv2.drawKeypoints(lena_rot45, kp2, lena_rot45, color=(255, 0, 255))

sift_cat1 = np.hstack((sift_original_lena, sift_lena_rot45))        # 對(duì)提取特征點(diǎn)后的圖像進(jìn)行橫向拼接
cv2.imwrite("sift_cat1.png", sift_cat1)
print('原圖與旋轉(zhuǎn)圖 特征點(diǎn)繪制圖像已保存')
cv2.imshow("sift_point1", sift_cat1)
cv2.waitKey()

# 特征點(diǎn)匹配
# K近鄰算法求取在空間中距離最近的K個(gè)數(shù)據(jù)點(diǎn)，并將這些數(shù)據(jù)點(diǎn)歸為一類
start = time.time()     # 計(jì)算匹配點(diǎn)匹配時(shí)間
bf = cv2.BFMatcher()
matches1 = bf.knnMatch(des1, des2, k=2)
print('用于 原圖和旋轉(zhuǎn)圖 圖像匹配的所有特征點(diǎn)數(shù)目：', len(matches1))

# 調(diào)整ratio
# ratio=0.4：對(duì)于準(zhǔn)確度要求高的匹配；
# ratio=0.6：對(duì)于匹配點(diǎn)數(shù)目要求比較多的匹配；
# ratio=0.5：一般情況下。
ratio1 = 0.5
good1 = []

for m1, n1 in matches1:
    # 如果最接近和次接近的比值大于一個(gè)既定的值，那么我們保留這個(gè)最接近的值，認(rèn)為它和其匹配的點(diǎn)為good_match
    if m1.distance < ratio1 * n1.distance:
        good1.append([m1])

end = time.time()
print("匹配點(diǎn)匹配運(yùn)行時(shí)間:%.4f秒" % (end-start))

# 通過對(duì)good值進(jìn)行索引，可以指定固定數(shù)目的特征點(diǎn)進(jìn)行匹配，如good[:20]表示對(duì)前20個(gè)特征點(diǎn)進(jìn)行匹配
match_result1 = cv2.drawMatchesKnn(original_lena, kp1, lena_rot45, kp2, good1, None, flags=2)
cv2.imwrite("match_result1.png", match_result1)

print('原圖與旋轉(zhuǎn)圖 特征點(diǎn)匹配圖像已保存')
print("=========================================")
print("=========================================")
print("原圖與旋轉(zhuǎn)圖匹配對(duì)的數(shù)目:", len(good1))

for i in range(2):
    print("匹配", i)
    print("數(shù)據(jù)類型:", type(good1[i][0]))
    print("描述符之間的距離:", good1[i][0].distance)
    print("查詢圖像中描述符的索引:", good1[i][0].queryIdx)
    print("目標(biāo)圖像中描述符的索引:", good1[i][0].trainIdx)

print("=========================================")
print("=========================================")
cv2.imshow("original_lena and lena_rot45 feature matching result", match_result1)
cv2.waitKey()

# 將待匹配圖像通過旋轉(zhuǎn)、變換等方式將其與目標(biāo)圖像對(duì)齊，這里使用單應(yīng)性矩陣。
# 單應(yīng)性矩陣有八個(gè)參數(shù)，如果要解這八個(gè)參數(shù)的話，需要八個(gè)方程，由于每一個(gè)對(duì)應(yīng)的像素點(diǎn)可以產(chǎn)生2個(gè)方程(x一個(gè)，y一個(gè))，那么總共只需要四個(gè)像素點(diǎn)就能解出這個(gè)單應(yīng)性矩陣。
if len(good1) > 4:
    ptsA = np.float32([kp1[m[0].queryIdx].pt for m in good1]).reshape(-1, 1, 2)
    ptsB = np.float32([kp2[m[0].trainIdx].pt for m in good1]).reshape(-1, 1, 2)
    ransacReprojThreshold = 4
    # RANSAC算法選擇其中最優(yōu)的四個(gè)點(diǎn)
    H, status =cv2.findHomography(ptsA, ptsB, cv2.RANSAC, ransacReprojThreshold)
    imgout = cv2.warpPerspective(lena_rot45, H, (original_lena.shape[1], original_lena.shape[0]),
                                 flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP)

    cv2.imwrite("imgout.png", imgout)
    cv2.imshow("lena_rot45's result after transformation", imgout)
    cv2.waitKey()