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

旋轉傻烏龜，效果就是將一隻烏龜在視窗中同時進行平移、縮放和旋轉，由於最後看起來樣子比較傻，因此得名“旋轉傻烏龜”。

效果視頻：

一、幾何變換的矩陣表示

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節的旋轉，可以表示為下麵形式。除此之外，還可以在此基礎上進行旋轉平移，只要在左邊依次乘上相應矩陣即可。

二、旋轉傻烏龜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)

三、旋轉傻烏龜自實現

 這個自己用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)