旋轉傻烏龜——幾何變換實踐

時間 2020-01-29

標籤旋轉烏龜幾何變換實踐简体版

原文原文鏈接

這兩天新型肺炎病例是指數上升啊！呆在家裏沒事幹，正好想起以前FPGA大賽上有個老哥作了一個圖像旋轉做品，還在羣裏發了技術報告。無聊之下就打算學習一下，而後就順便把平移、旋轉、縮放這些幾何變換都看了，最後決定把這三個綜合起來寫個「旋轉傻烏龜」的動畫。先是用OpenCV內置函數實現了下，感受不過癮，又本身寫了一遍。老規矩，仍是把學過的、作過的東西記錄下來！html

旋轉傻烏龜，效果就是將一隻烏龜在窗口中同時進行平移、縮放和旋轉，因爲最後看起來樣子比較傻，所以得名「旋轉傻烏龜」。算法

效果視頻：編程

1、幾何變換的矩陣表示

1.1　平移的表示

上圖中的三種表示方法第二種是OpenCV要求的方式，但第一種形式表示起來更具統一性，所以我更傾向於第一種。但不管哪種，都能展開成第三種的形式。第三種很是直觀的反映了平移，只是須要注意正負號的選取——在編程中，圖像通常以左上角爲（0，0）點。這也就是說，創建座標系的時候，X軸以右正方向，Y軸如下爲正方向。以上矩陣表示將圖像向右平移x0，向下平移y0，也能夠認爲是將座標系向左平移x0，向上平移y0。平移能夠形象地表示以下：編程語言

1.2　以左上角爲定點縮放的表示

縮放最容易理解，就是將橫縱座標乘以縮放比例。因爲咱們以左上角爲座標系原點，因此左上角點的位置並不會變化。函數

1.3　以左上角點爲中心旋轉的表示

在本文中，規定順時針方向旋轉，θ爲正；逆時針旋轉，θ爲負。旋轉先後的座標關係推導也不難，以下圖所示，旋轉前先求出旋轉半徑L，旋轉後根據L求出座標。性能

爲了以後表述的簡潔，咱們將這三節中的矩陣分別用特定符號簡記：學習

1.4　以任一點爲中心旋轉的表示

有了以上的基礎，咱們就能夠研究更加複雜的變換。例如咱們想以任一點(x0,y0)爲中心旋轉，而咱們推導的R(θ)只適用於以座標系原點爲中心旋轉。所以，咱們能夠將圖像向上平移x0，向左平移y0，使(x0,y0)點平移到座標系原點；而後再旋轉，旋轉完後再向下平移x0，向右平移y0回到原來位置，這一過程可用三個基礎基礎矩陣表示成以下形式，注意三個矩陣順序不能調換。測試

1.5　以任一點爲定點縮放的表示

方法同1.4節的旋轉，能夠表示爲下面形式。除此以外，還能夠在此基礎上進行旋轉平移，只要在左邊依次乘上相應矩陣便可。優化

2、旋轉傻烏龜OpenCV函數實現

OpenCV提供了仿射變換函數warpAffine。在輸入參數中，M表示變換矩陣，能夠是平移、旋轉和縮放矩陣等；dsize是輸入圖像的大小；flags是插值方式，通常採用默認的雙線性插值。動畫

至於M的獲取，平移矩陣只能本身構造；二旋轉矩陣能夠由函數getRotationMatrix2D獲得。輸入參數中，center表示旋轉中心的座標；angle爲旋轉角度，逆時針爲正；scale是縮放比例。可見這個函數同時包攬了旋轉和縮放的功能。

個人思路是，用正弦函數生成一系列軌跡點，烏龜每到達一個軌跡點，就旋轉必定角度，縮放必定比例，而軌跡點的跟蹤就是烏龜中心的平移。根據以前的說的原理，咱們先讓整個圖像繞自身中心旋轉和縮放，縮放後的烏龜應該是在整個圖像的中間，爲了讓它中心和軌跡重合，就使用平移變換，此時平移的距離應該是path-center。整個過程的代碼以下：

 1 import cv2
 2 import numpy as np
 3 import time
 4 
 5 img = cv2.imread('image/turtle.jpg')
 6 size = img.shape[:-1]
 7 cv2.namedWindow('img')
 8 
 9 #平移矩陣
10 def GetMoveMatrix(x,y):
11     M = np.zeros((2, 3), dtype=np.float32)
12 
13     M.itemset((0, 0), 1)
14     M.itemset((1, 1), 1)
15     M.itemset((0, 2), x)
16     M.itemset((1, 2), y)
17 
18     return M
19 
20 if __name__ == '__main__':
21 
22     # shape和座標是顛倒的
23     center_x = size[1]/2
24     center_y = size[0]/2
25     #計時
26     start_time = time.time()
27 
28     for x in np.linspace(0,2*np.pi,100):
29         #角度、縮放
30         angle = -360*x/2/np.pi
31         scale = 0.2+0.2*np.sin(x)
32         #軌跡
33         path_x = x*50+100
34         path_y = (np.sin(x)+1)*100+100
35         #旋轉、平移矩陣
36         M1 = cv2.getRotationMatrix2D((center_x, center_y), angle, scale)
37         M2 = GetMoveMatrix(path_x-center_x,path_y-center_y)
38         #仿射變換
39         rotate = cv2.warpAffine(img,M1,size)
40         dst = cv2.warpAffine(rotate,M2,size)
41 
42         # cv2.imshow('img',dst)
43         # cv2.waitKey(1)
44     #花費125ms
45     print(time.time()-start_time)

3、旋轉傻烏龜自實現

這個本身用Python實現的話，性能就至關重要了，尤爲是雙線性插值，若是不優化的話，慢得簡直可讓你懷疑人生。好比，通常的是用兩個for循環迭代，代碼以下。在這個項目裏，這個函數執行一次須要花費1.4s的時間。因此不優化的話，這隻烏龜真的是名副其實了！

 1 def InterLinearMap(img,size,mapx,mapy):
 2 
 3     dst = np.zeros(img.shape,dtype=np.uint8)
 4 
 5     for row in range(size[0]):
 6         for col in range(size[1]):
 7 
 8             intx = np.int32(mapx.item(row,col))
 9             inty = np.int32(mapy.item(row,col))
10             partx = mapx.item(row,col)-intx
11             party = mapy.item(row,col)-inty
12             resx = 1-partx
13             resy = 1-party
14 
15             if party==0 and partx==0:
16                 result=img[inty,intx]
17             else:
18                 result = ((img[inty,intx]*resx+img[inty,intx+1]*partx)*resy
19                           +(img[inty+1,intx]*resx+img[inty+1,intx+1]*partx)*party)
20 
21             dst[row,col]=np.uint8(result+0.5)
22 
23     return dst

那怎麼辦？網上有一些優化的方法，主要是將浮點運算轉成整數運算，這個方法對於FPGA這樣的邏輯器件最適合不過了——但別忘了，我如今用的是Python，整數運算實際上也會被轉成浮點運算，因此這個方法顯然不適用。我採用的優化是進行矩陣化，據我所知，不少編程語言只要是支持矩陣運算的，其運算都是優化過的。對於雙線性插值和仿射變換，運用矩陣也是很合適，只是寫起來會有點抽象。。。

首先，先把生成變換矩陣的函數寫出來，代碼以下。要注意numpy的三角函數接受的參數是弧度制。

 1 #縮放矩陣
 2 def GetResizeMatrix(scalex,scaley):
 3     M = np.zeros((3,3),dtype=np.float32)
 4 
 5     M.itemset((0,0),scalex)
 6     M.itemset((1,1),scaley)
 7     M.itemset((2,2),1)
 8 
 9     return M
10 #平移矩陣
11 def GetMoveMatrix(x,y):
12     M = np.zeros((3, 3), dtype=np.float32)
13 
14     M.itemset((0, 0), 1)
15     M.itemset((1, 1), 1)
16     M.itemset((2, 2), 1)
17     M.itemset((0, 2), x)
18     M.itemset((1, 2), y)
19 
20     return M
21 #旋轉矩陣
22 def GetRotationMatrix(angle):
23     M = np.zeros((3, 3), dtype=np.float32)
24 
25     M.itemset((0, 0), np.cos(angle))
26     M.itemset((0, 1), -np.sin(angle))
27     M.itemset((1, 0), np.sin(angle))
28     M.itemset((1, 1), np.cos(angle))
29     M.itemset((2, 2), 1)
30 
31     return M

接下來寫仿射變換函數，輸入參數爲圖片數據、變換矩陣和輸入圖片的大小。這裏應該要有逆向思惟——如今我要獲得變換後的圖片，就是要求各座標位置上的色彩，而色彩取樣自變換前圖像上的一點（這點的座標可能不是整數），也就是說咱們要將變換後的座標映射到變換前的座標。再來看以前的公式（下圖左，爲了方便，將變換矩陣合成爲一個矩陣A），如今咱們已知的是左邊部分，而要求的映射是等式右邊的XY，所以咱們將A拿到左邊，獲得另外一個公式（下圖右），並依據這個公式，寫出仿射變換函數。

 1 def WarpAffine(img,Mat,size):
 2 
 3     rows = size[0]
 4     cols = size[1]
 5     #生成矩陣[X Y 1]
 6     ones = np.ones((rows, cols), dtype=np.float32)
 7     #gridx/gridy -> shape(rows,cols)
 8     gridx,gridy= np.meshgrid(np.arange(0, cols),np.arange(0, rows))
 9     #dst -> shape(3,rows,cols)
10     dst = np.stack((gridx, gridy, ones))
11 
12     #求逆矩陣 M -> shape(3,3)
13     Mat = np.linalg.inv(Mat)
14     #得到矩陣[x,y,1] -> shape(3,rows,cols)
15     src = np.tensordot(Mat,dst,axes=[[-1],[0]])
16     
17     #mapx/mapy -> shape(rows,cols)
18     mapx = src[0]#座標非整數
19     mapy = src[1]#座標非整數
20     #仿射出界的設爲原點
21     flags = (mapy > rows - 2) + (mapy < 0) + (mapx > cols - 2) + (mapx < 0)
22     mapy[flags] = 0
23     mapx[flags] = 0
24     #雙線性插值
25 
26     result = InterLinearMap(img, mapx, mapy)
27 
28     return result

再解決雙線性插值，關於該算法的原理挺簡單的，讀者能夠網上查找（提一點，理解雙線性插值時能夠想象3D模型，Z軸爲灰度值）。對於該函數，借鑑一下remap函數，輸入參數設兩個map，分別表示x,y的映射。map的大小跟圖片大小相同，也就是說，一共有rows*cols點須要插值，除了用兩個for迭代，咱們也能夠將rows和cols做爲矩陣的兩個額外維度，表示樣本數。計算的話，利用矩陣的點乘代替凌亂的長算式，顯得很簡潔，公式以下：

代碼以下，經測試，執行一次該函數，花費時間爲45ms，這要比原來的1.4s快多了（實在不知道該怎麼進一步優化了，mxy、img下表索引、求和各花了15ms）

def InterLinearMap(img,mapx,mapy):

    #(rows,cols)
    inty = np.int32(mapy)
    intx = np.int32(mapx)
    nxty = 1+inty
    nxtx = 1+intx
    #(rows,cols)
    party = mapy - inty
    partx = mapx - intx
    resy = 1-party
    resx = 1-partx

    #(4,rows,cols)
    mxy = np.stack((resy*partx,resy*resx,partx*party, resx*party))
    mxy = np.expand_dims(mxy,axis=-1)

    #(4,rows,cols,3)
    mf = np.stack((img[inty,nxtx],img[inty,intx],img[nxty,nxtx],img[nxty,intx]))

    #res -> shape(rows,cols,3)
    res = np.sum(mxy*mf,axis=0)
    res = np.uint8(res+0.5)

    return res

綜上，給出完整代碼：

import cv2
import numpy as np

img = cv2.imread('image/turtle.jpg')
size = img.shape[:-1]
cv2.namedWindow('img')

#縮放矩陣
def GetResizeMatrix(scalex,scaley):
    M = np.zeros((3,3),dtype=np.float32)

    M.itemset((0,0),scalex)
    M.itemset((1,1),scaley)
    M.itemset((2,2),1)

    return M
#平移矩陣
def GetMoveMatrix(x,y):
    M = np.zeros((3, 3), dtype=np.float32)

    M.itemset((0, 0), 1)
    M.itemset((1, 1), 1)
    M.itemset((2, 2), 1)
    M.itemset((0, 2), x)
    M.itemset((1, 2), y)

    return M
#旋轉矩陣
def GetRotationMatrix(angle):
    M = np.zeros((3, 3), dtype=np.float32)

    M.itemset((0, 0), np.cos(angle))
    M.itemset((0, 1), -np.sin(angle))
    M.itemset((1, 0), np.sin(angle))
    M.itemset((1, 1), np.cos(angle))
    M.itemset((2, 2), 1)

    return M

def InterLinearMap(img,mapx,mapy):

    #(rows,cols)
    inty = np.int32(mapy)
    intx = np.int32(mapx)
    nxty = 1+inty
    nxtx = 1+intx
    #(rows,cols)
    party = mapy - inty
    partx = mapx - intx
    resy = 1-party
    resx = 1-partx

    #(4,rows,cols)
    mxy = np.stack((resy*partx,resy*resx,partx*party, resx*party))
    mxy = np.expand_dims(mxy,axis=-1)

    #(4,rows,cols,3)
    mf = np.stack((img[inty,nxtx],img[inty,intx],img[nxty,nxtx],img[nxty,intx]))

    #res -> shape(rows,cols,3)
    res = np.sum(mxy*mf,axis=0)
    res = np.uint8(res+0.5)

    return res



def WarpAffine(img,Mat,size):

    rows = size[0]
    cols = size[1]
    #生成矩陣[X Y 1]
    ones = np.ones((rows, cols), dtype=np.float32)
    #gridx/gridy -> shape(rows,cols)
    gridx,gridy= np.meshgrid(np.arange(0, cols),np.arange(0, rows))
    #dst -> shape(3,rows,cols)
    dst = np.stack((gridx, gridy, ones))

    #求逆矩陣 M -> shape(3,3)
    Mat = np.linalg.inv(Mat)
    #得到矩陣[x,y,1] -> shape(3,rows,cols)
    src = np.tensordot(Mat,dst,axes=[[-1],[0]])
    
    #mapx/mapy -> shape(rows,cols)
    mapx = src[0]#座標非整數
    mapy = src[1]#座標非整數
    #仿射出界的設爲原點
    flags = (mapy > rows - 2) + (mapy < 0) + (mapx > cols - 2) + (mapx < 0)
    mapy[flags] = 0
    mapx[flags] = 0
    #雙線性插值

    result = InterLinearMap(img, mapx, mapy)

    return result



if __name__ == '__main__':

    center_x = size[1]/2
    center_y = size[0]/2

    for x in np.linspace(0,2*np.pi,100):

        angle = 360*x/2/np.pi
        scale = 0.2+0.2*np.sin(x)

        path_x = x*50+100
        path_y = (np.sin(x)+1)*100+100

        M = GetMoveMatrix(path_x,path_y)@GetRotationMatrix(x)\
            @GetResizeMatrix(scale,scale)@GetMoveMatrix(-center_x,-center_y)

        dst = WarpAffine(img,M,size)
        cv2.imshow('img',dst)
        cv2.waitKey(1)