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