計算機視覺2D幾何基元及其變換介紹和OpenCV WarpPerspective源碼分析

2D圖像幾何基元

---

通常的，表示一個2d幾何基元只用兩個維度(好比x，y)就能夠表示了，可是在計算機視覺研究中，爲了統一對2d幾何基元的操做（後面講到的仿射，透射變換）,通常會以增廣矢量的方式表示幾何基元。

齊次座標將本來n維的座標用一個n+1維的座標表示，其兩個基本做用爲：

    1. 區分n維空間的點和向量，一個點的第n+1維爲非零值，而向量的n+1維爲0

    2. 統一幾何基元的旋轉，平移，拉伸，投影等操做（只用一個矩陣就能夠表示）

2D點：2D點的齊次座標表示爲：，其中僅在尺度上不一樣的矢量被視爲等同的，被稱做爲2D投影空間，其次矢量能夠經過除以最後你一個元素來轉換爲非齊次矢量X，即：其中，是增廣矢量，若是最後w爲0，則稱此點爲理想點或者無窮遠點，它沒有等同的非齊次表達

2D直線：2D直線能夠用齊次向量表示，其對應的直線方程爲：，其規範化表達爲，若是 l = (0,0,1)，則包含全部無窮遠點

2D圓錐曲線：

 

2D圖像變換及示例

---

利用齊次座標，咱們能夠把統一2D變換操做

平移變換：

，其中I是2x2的單位矩陣

旋轉變換：

放縮變換（x，y分別放縮）：

仿射變換：

投影變換：

由於而仿射變換是投影變化的子集，全部對2D幾何基元的操做都以用一個3x3的Homography表示，這樣一來對圖像的屢次2D變換就至關於圖像對多個Homography矩陣的乘積，即R = H1*H2*H3...Hn*r 

 

咱們用Python代碼表示這個過程

 1 import numpy as np
 2 import cv2 as cv
 3 from numba import jit
 4 import matplotlib.pyplot as plt
 5 from matplotlib.font_manager import FontProperties
 6 
 7 def MoveT(tx, ty, lastOp = None)->np.ndarray:
 8     op = np.array([[1, 0, tx], [0, 1, ty], [0, 0, 1]], dtype=np.float32)
 9     if lastOp is not None:
10         op = np.dot(op, lastOp)
11     return op
12 
13 
14 def RotateT(r, lastOp = None)->np.ndarray:
15     op = np.array([[np.cos(r), -np.sin(r), 0], [np.sin(r), np.cos(r), 0], [0, 0, 1]], dtype=np.float32)
16     if lastOp is not None:
17         op = np.dot(op, lastOp)
18     return op
19 
20 
21 def ZoomT(rx, ry, lastOp = None)->np.ndarray:
22     op = np.array([[rx, 0, 0], [0, ry, 0], [0, 0, 1]], dtype=np.float32)
23     if lastOp is not None:
24         op = np.dot(op, lastOp)
25     return op
26 
27 def Transform(imgSrc:np.ndarray, op)->np.ndarray:
28     h, w = imgSrc.shape[:2]
29     imgDst = np.zeros(imgSrc.shape, dtype=imgSrc.dtype)
30     for i in range(h):
31         v1 = np.stack((np.arange(w), np.ones(w)*i, np.ones(w)),axis=-1)
32         v2 = np.dot(v1, op)
33         tpx, tpy, tpz = np.hsplit(v2, 3)
34         for iy, ix, iz, j in zip(tpy, tpx, tpz, range(w)):
35             py, px = int(iy/iz), int(ix/iz)
36             if 0<= py < h and 0 <= px < w:
37                 imgDst[int(py), int(px)] = imgSrc[i, j]
38 
39     return imgDst

1 if __name__ == "__main__":
2     font_set = FontProperties(fname=r"c:\windows\fonts\msyh.ttc", size=12)
3     imgSrc = plt.imread("E:/Users/Administrator/pictures/Test/user.jpg")
4     op = np.transpose(MoveT(10,30, RotateT(np.pi/12, ZoomT(1.1, 1.2))))
5 
6     imgDst = Transform(imgSrc, op)
7     plt.figure(1), plt.imshow(imgDst), plt.title("Resuult", fontproperties=font_set)
8     plt.show()

 

咱們對圖進行了以下操做：

1. 對y尺度放大1.1倍，x尺度放大1.2倍

2. 繞z軸（垂直於x,y）旋轉π/12度

3. x方向移動30像素，y方向移動10像素

fig1: 原圖	fig2: 仿射變換結果

 

能夠看到這個結果並非很好，在結果中發現一些很明顯的空洞，由於咱們對圖像進行了放大操做和旋轉操做，在這一個過程當中結果圖像的有一些像素的值並不能從仿射變換中得出，從而造成空洞，可見直接直接對圖像進行正向變換效果是很差的。下面咱們將會看到利用反向變換+插值的方法解決空洞的問題。

 

 

利用透射變換來進行傾斜校訂

---

透射變換的一個應用就是拿來校訂圖像，原理很簡單，就是把根據關鍵點找到變換矩陣

 

假如使用矩陣，而且設h33 = 1，那麼就要把剩下的8個參數給解出來，8個未知數須要8個方程，又由於從仿射變換的對應關係，能夠列方程爲：

也就是說只要找8個點的以及其兩兩對應關係，就能夠找到透射變換的矩陣，也就能夠對圖像進行校訂了

透射變換有放縮和旋轉等操做，正如上一小節所說，若是僅使用正向變換，那麼確定會有空洞現象的發生。若是使用反向變換，那麼這個問題也就迎刃而解了，咱們要作的就是掃描反向變換的位置，在這裏咱們能夠採用原圖像正向變換後上下左右的邊界點肯定掃描區域，而後一個一個點反向變換回原圖像中，並根據原圖像來進行插值

 

Python代碼

  1 # WarpCorrection.py
  2 
  3 @jit
  4 def BilinearInterpolation(imgSrc:np.ndarray, h, w, sx:float, sy:float)->float:
  5     """
  6     對圖片的指定位置作雙線性插值
  7     :param imgSrc:源圖像
  8     :param h: src的高度
  9     :param w: src的寬度
 10     :param sx: x位置
 11     :param sy: y位置
 12     :return: 所插入的值
 13     """
 14     intSx, intSy = int(sx), int(sy)
 15     if 0 <= intSx  < w - 1 and 0 <= intSy < h - 1:
 16         x1, x2 = intSx, intSx + 1
 17         y1, y2 = intSy, intSy + 1
 18         H1 = np.dot(np.array([x2 - sx, sx - x1]), imgSrc[y1: y2 + 1, x1:x2 + 1])
 19         return H1[0]*(y2 - sy) + H1[1]*(sy - y1)
 20     else:
 21         return imgSrc[intSy, intSx]
 22 
 23 def WarpCorrection(imgSrc:np.ndarray, dots:tuple)->np.ndarray:
 24     assert len(dots) == 4
 25 
 26     # 四個點的順序必定要按照左上，右上，右下，左下的順時針順序點
 27     d1, d2, d3, d4 = dots
 28     x1, x2, x3, x4 = d1[0], d2[0], d3[0], d4[0]
 29     y1, y2, y3, y4 = d1[1], d2[1], d3[1], d4[1]
 30     assert x1 < x2
 31     assert x4 < x3
 32     assert y1 < y4
 33     assert y2 < y3
 34 
 35     objW = np.round(np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2))
 36     objH = np.round(np.sqrt((x2 - x3) ** 2 + (y2 - y3) ** 2))
 37 
 38     # 在這裏我簡單地設爲把所輸入的四個點的位置，經過2D變換，變換爲長方形的四個頂點的位置（以x1爲起點）
 39     t1, t2, t3, t4 = (y1, x1), (y1, x1 + objW), (y1 + objH, x1 + objW), (y1 + objH, x1),
 40 
 41     rx1, rx2, rx3, rx4 = t1[1], t2[1], t3[1], t4[1]
 42     ry1, ry2, ry3, ry4 = t1[0], t2[0], t3[0], t4[0]
 43 
 44     # ================Step 0: 根據 8個點兩兩對應關係找到Homography矩陣================
 45     # 把8個約束寫成方程組，以矩陣的形式表達
 46     m = np.array([
 47                   [y1, x1, 1, 0, 0, 0, -ry1 * y1, -ry1 * x1],
 48                   [0, 0, 0, y1, x1, 1, -rx1 * y1, -rx1 * x1],
 49                   [y2, x2, 1, 0, 0, 0, -ry2 * y2, -ry2 * x2],
 50                   [0, 0, 0, y2, x2, 1, -rx2 * y2, -rx2 * x2],
 51                   [y3, x3, 1, 0, 0, 0, -ry3 * y3, -ry3 * x3],
 52                   [0, 0, 0, y3, x3, 1, -rx3 * y3, -rx3 * x3],
 53                   [y4, x4, 1, 0, 0, 0, -ry4 * y4, -ry4 * x4],
 54                   [0, 0, 0, y4, x4, 1, -rx4 * y4, -rx4 * x4],
 55                 ])
 56 
 57     vectorSrc = np.array([ry1, rx1, ry2, rx2, ry3, rx3, ry4, rx4])
 58     vectorSrc.shape = (1, 8)
 59     HFlat = np.dot(np.linalg.inv(m), np.transpose(vectorSrc))
 60     a, b, c, d, e, f, g, h = HFlat[0, 0],HFlat[1, 0],HFlat[2, 0],HFlat[3, 0],HFlat[4, 0],HFlat[5, 0],HFlat[6, 0],HFlat[7, 0]
 61 
 62     H = np.array([[a, b, c],
 63                   [d, e, f],
 64                   [g, h, 1]], dtype=np.float32)
 65 
 66     # ================Step 1: 經過對原圖像四個頂點進行正向投射變換，肯定目標圖像區域================
 67     height, width = imgSrc.shape[:2]
 68     matrixOriginVertex = np.array([[0, 0, 1],
 69                                    [0, width - 1, 1],
 70                                    [height - 1, width - 1, 1] ,
 71                                    [height - 1, 0, 1]])
 72 
 73     result = np.dot(matrixOriginVertex, np.transpose(H))
 74     minX = int(min(result[0, 1]/result[0, 2], result[1, 1]/result[1, 2], result[2, 1]/result[2, 2], result[3, 1]/result[3, 2]))
 75     maxX = int(max(result[0, 1]/result[0, 2], result[1, 1]/result[1, 2], result[2, 1]/result[2, 2], result[3, 1]/result[3, 2]))
 76     minY = int(min(result[0, 0]/result[0, 2], result[1, 0]/result[1, 2], result[2, 0]/result[2, 2], result[3, 0]/result[3, 2]))
 77     maxY = int(max(result[0, 0]/result[0, 2], result[1, 0]/result[1, 2], result[2, 0]/result[2, 2], result[3, 0]/result[3, 2]))
 78 
 79     # ================Step 2: 反向變換+雙二次插值校訂圖像================
 80     vtr = np.empty((0,3),dtype=np.float32)
 81     for i in range(minY, maxY):
 82         arr1 = np.arange(minX, maxX)
 83         arr2 = np.ones(maxX - minX)
 84         vt1 = np.stack((arr2*i, arr1 , arr2), axis=-1)
 85         vtr = np.concatenate((vtr, vt1), axis=0)
 86 
 87     # 請注意，由於傳進去的是規範化後(Y, X, 1)的值，因此獲得的實際上是(y/Z, x/Z, 1/Z的值)
 88     vts = np.dot(vtr,np.linalg.inv(np.transpose(H)))
 89     dstHeight, dstWidth = maxY - minY + 1, maxX - minX + 1
 90     imgDst = np.zeros((dstHeight, dstWidth, imgSrc.shape[2]), dtype=imgSrc.dtype)
 91 
 92     for (r, s) in zip(vtr, vts):
 93         ry, rx = int(r[0]), int(r[1])
 94         iy, ix = s[:2]
 95         # 須要解 [y, x] = [iy*(g*y + h*x + 1), ix*(g*y + h*x + 1)]這個方程
 96         TH = np.linalg.inv(np.array([[iy * g - 1, iy * h],
 97                                      [ix * g, ix * h - 1]]))
 98 
 99         vxy = np.dot(TH, np.array([[-iy], [-ix]]))
100         sy, sx = vxy[0, 0], vxy[1, 0]
101 
102         if 0 <= round(sy) < height and 0 <= round(sx) < width:
103             imgDst[ry - minY, rx - minX] = BilinearInterpolation(imgSrc, height, width, sx, sy)
104 
105     return imgDst

這裏使用了雙二次插值，雙二次插值的公式推導以下：

fig3. 雙二次插值公式推導	fig4. 雙二次插值圖像解釋

上面右圖所展現的那樣，雙二次插值其實就是一個加權平均操做，爲了和雙二次插值的效果作對比，我加了個最鄰近插值

 

繼續往程序裏面加點東西，來測試下實際效果，在這裏我使用了OpenCV的Highgui的包

 1 # -*- coding: utf-8 -*-
 2 import matplotlib.pyplot as plt
 3 from matplotlib.font_manager import FontProperties
 4 from WarpTransform import *
 5 from multiprocessing.dummy import Process
 6 
 7 windowCount = int(0)
 8 mainWinName = "source"
 9 
10 def WarpImage(imgSrc:np.ndarray, dots:tuple, count)->None:
11     imgBLiner, imgNearest = WarpCorrection(imgSrc, dots)
12 
13     winName:str = f"result BLiner {count}"
14     cv.namedWindow(winName)
15     cv.imshow(winName, imgBLiner)
16 
17     winName:str = f"result nearest {count}"
18     cv.namedWindow(winName)
19     cv.imshow(winName, imgNearest)
20 
21     cv.waitKey(0)
22     cv.destroyWindow(winName)
23 
24 class WarpCorrectionMgr:
25     def __init__(self, imgSrc):
26         self.__clickTime = 0
27         self.__imgSrc = imgSrc.copy()
28         self.__imgDrawn = imgSrc.copy()
29         self.__dots = []
30 
31     @property
32     def sourceImage(self):
33         return self.__imgSrc
34 
35     @property
36     def drawnImage(self):
37         return self.__imgDrawn
38 
39     @drawnImage.setter
40     def drawnImage(self, newImg):
41         self.__imgDrawn = newImg
42 
43     @property
44     def clickTime(self):
45         return self.__clickTime
46 
47     @clickTime.setter
48     def clickTime(self, v):
49         self.__clickTime = v
50 
51     @property
52     def dots(self):
53         return self.__dots
54 
55     @staticmethod
56     def MouseCallback(event, x, y, flags, param):
57         # 四個點的順序必定要按照左上，右上，右下，左下的順時針順序點
58         if event == cv.EVENT_LBUTTONDBLCLK:
59             clickTime = param.clickTime
60             cv.circle(param.drawnImage, (x, y), 8, (0, 0, 255),-1)
61             param.dots.append((x, y))
62             cv.imshow(mainWinName, param.drawnImage)
63 
64             if clickTime + 1 == 4:
65                 global windowCount
66                 p = Process(target=WarpImage, args=(param.sourceImage, param.dots.copy(), windowCount))
67                 p.daemon = True
68                 p.start()
69 
70                 param.drawnImage = param.sourceImage.copy()
71                 cv.imshow(mainWinName,param.sourceImage)
72                 param.dots.clear()
73                 windowCount += 1
74 
75             param.clickTime = (clickTime + 1) % 4
76 
77 
78 if __name__ == "__main__":
79     cv.namedWindow(mainWinName)
80     imgSrc = cv.imread("E:/Users/Administrator/pictures/Test/skew.jpg")
81     imgSrc = cv.resize(imgSrc, (int(imgSrc.shape[1]/4), int(imgSrc.shape[0]/4)))
82 
83     mgr = WarpCorrectionMgr(imgSrc)
84     cv.setMouseCallback(mainWinName, WarpCorrectionMgr.MouseCallback, mgr)
85 
86     cv.imshow(mainWinName, imgSrc)
87     cv.waitKey(0)
88     cv.destroyAllWindows()

爲了對比最近鄰插值和雙二次插值的差異，我把WarpCorrection改了一下：

    # 請注意，由於傳進去的是規範化後(Y, X, 1)的值，因此獲得的實際上是(y/Z, x/Z, 1/Z的值)
    vts = np.dot(vtr,np.linalg.inv(np.transpose(H)))
    dstHeight, dstWidth = maxY - minY + 1, maxX - minX + 1
    imgBLiner = np.zeros((dstHeight, dstWidth, imgSrc.shape[2]), dtype=imgSrc.dtype)
    imgNearest = np.zeros((dstHeight, dstWidth, imgSrc.shape[2]), dtype=imgSrc.dtype)
    
    for (r, s) in zip(vtr, vts):
        ry, rx = int(r[0]), int(r[1])
        iy, ix = s[:2]
        # 須要解 [y, x] = [iy*(g*y + h*x + 1), ix*(g*y + h*x + 1)]這個方程
        TH = np.linalg.inv(np.array([[iy * g - 1, iy * h],
                                     [ix * g, ix * h - 1]]))

        vxy = np.dot(TH, np.array([[-iy], [-ix]]))
        sy, sx = vxy[0, 0], vxy[1, 0]

        if 0 <= round(sy) < height and 0 <= round(sx) < width:
            imgBLiner[ry - minY, rx - minX] = BilinearInterpolation(imgSrc, height, width, sx, sy)
            imgNearest[ry - minY, rx - minX] = imgSrc[int(round(sy)),int(round(sx))]

    return imgBLiner, imgNearest

運行下程序來看下效果，上面程序的做用就是選中要校訂圖像的四個角，而後校訂（PS：實際運行可能會至關慢，我沒作優化）

原圖

Fig3. 最近鄰插值（邊緣很粗糙，看「清華大學出版社」幾個字）	Fig4. 雙二次插值（毛刺不多，可是圖像模糊了）

OpenCV庫中WarpPerspective源代碼

---

其實OpenCV是帶了仿射變換和透射變換的API的，先來看下投影變換的源碼（在imgwarp.cpp）（部分關鍵代碼，Assert和IPP優化這裏不展現了）：

 

 1 void cv::warpPerspective( InputArray _src, OutputArray _dst, InputArray _M0,
 2                           Size dsize, int flags, int borderType, const Scalar& borderValue )
 3 {
 4     Mat src = _src.getMat(), M0 = _M0.getMat();
 5     _dst.create( dsize.area() == 0 ? src.size() : dsize, src.type() );
 6     Mat dst = _dst.getMat();
 7 
 8     if( dst.data == src.data )
 9         src = src.clone();
10 
11     double M[9];
12     Mat matM(3, 3, CV_64F, M);
13     int interpolation = flags & INTER_MAX;
14     
15     // 插入方法不支持INTER_AREA，用INTER_LINEAR代替
16     if( interpolation == INTER_AREA )
17         interpolation = INTER_LINEAR;
18 
19     CV_Assert( (M0.type() == CV_32F || M0.type() == CV_64F) && M0.rows == 3 && M0.cols == 3 );
20     M0.convertTo(matM, matM.type());
21     
22     //這裏源碼是這樣操做的，若是能夠執行IPP，那麼就直接執行IPP指令加塊操做
23     //若是沒有說明傳入的變換矩陣已經轉置過了，那麼執行一次矩陣轉置
24     if( !(flags & WARP_INVERSE_MAP) )
25         invert(matM, matM);
26 
27     hal::warpPerspectve(src.type(), src.data, src.step, src.cols, src.rows, dst.data, dst.step, dst.cols, dst.rows,
28                         matM.ptr<double>(), interpolation, borderType, borderValue.val);
29 }

 

 1 namespace hal {
 2 
 3 void warpPerspectve(int src_type,
 4                     const uchar * src_data, size_t src_step, int src_width, int src_height,
 5                     uchar * dst_data, size_t dst_step, int dst_width, int dst_height,
 6                     const double M[9], int interpolation, int borderType, const double borderValue[4])
 7 {
 8     Mat src(Size(src_width, src_height), src_type, const_cast<uchar*>(src_data), src_step);
 9     Mat dst(Size(dst_width, dst_height), src_type, dst_data, dst_step);
10 
11     Range range(0, dst.rows);
12     WarpPerspectiveInvoker invoker(src, dst, M, interpolation, borderType, Scalar(borderValue[0], borderValue[1], borderValue[2], borderValue[3]));
13     
14     //這裏是一個opencv執行並行處理的一個方法，真正執行處理的代碼在WarpPerspectiveInvoker的函數對象操做符中
15     parallel_for_(range, invoker, dst.total()/(double)(1<<16));
16 }
17 
18 } // hal::

 

  1 class WarpPerspectiveInvoker :
  2     public ParallelLoopBody
  3 {
  4 public:
  5     WarpPerspectiveInvoker(const Mat &_src, Mat &_dst, const double *_M, int _interpolation,
  6                            int _borderType, const Scalar &_borderValue) :
  7         ParallelLoopBody(), src(_src), dst(_dst), M(_M), interpolation(_interpolation),
  8         borderType(_borderType), borderValue(_borderValue)
  9     {
 10         //注意M就是轉置矩陣
 11     }
 12 
 13     virtual void operator() (const Range& range) const
 14     {
 15         const int BLOCK_SZ = 32;
 16         short XY[BLOCK_SZ*BLOCK_SZ*2], A[BLOCK_SZ*BLOCK_SZ];
 17         int x, y, x1, y1, width = dst.cols, height = dst.rows;
 18 
 19         int bh0 = std::min(BLOCK_SZ/2, height);            
 20         int bw0 = std::min(BLOCK_SZ*BLOCK_SZ/bh0, width);  //操做塊的寬
 21         bh0 = std::min(BLOCK_SZ*BLOCK_SZ/bw0, height);     //操做塊的高
 22 
 23         //下面的意思是，若是支持SSE4指令，那麼就先對轉置矩陣進行預處理，下面會繼續用到這個pwarp_impl_sse4
 24         #if CV_TRY_SSE4_1
 25         Ptr<opt_SSE4_1::WarpPerspectiveLine_SSE4> pwarp_impl_sse4;
 26         if(CV_CPU_HAS_SUPPORT_SSE4_1)
 27             pwarp_impl_sse4 = opt_SSE4_1::WarpPerspectiveLine_SSE4::getImpl(M);
 28         #endif
 29 
 30         //opencv對投射變換的操做是基於塊操做的，主要是爲了方即可以直接用SSE4指令給操做提速
 31         for( y = range.start; y < range.end; y += bh0 )
 32         {
 33             for( x = 0; x < width; x += bw0 )
 34             {
 35                 int bw = std::min( bw0, width - x);
 36                 int bh = std::min( bh0, range.end - y); // height
 37 
 38                 //_XY是一個雙通道的Mat，下面remap會用到
 39                 Mat _XY(bh, bw, CV_16SC2, XY), matA;
 40                 Mat dpart(dst, Rect(x, y, bw, bh));
 41 
 42                 for( y1 = 0; y1 < bh; y1++ )
 43                 {
 44                     short* xy = XY + y1*bw*2;
 45                     //就是[X,Y,W] = [x,y,1]*M^T的操做
 46                     //注意如今是分塊操做，因此纔有+ M[1]*(y + y1)這一項
 47                     double X0 = M[0]*x + M[1]*(y + y1) + M[2];
 48                     double Y0 = M[3]*x + M[4]*(y + y1) + M[5];
 49                     double W0 = M[6]*x + M[7]*(y + y1) + M[8];
 50 
 51                     //最近鄰插值
 52                     if( interpolation == INTER_NEAREST )
 53                     {
 54                         x1 = 0;
 55 
 56                         //這裏就是剛纔上面所說的pwarp_impl_sse4所用到的位置
 57                         #if CV_TRY_SSE4_1
 58                         if (pwarp_impl_sse4)
 59                             pwarp_impl_sse4->processNN(M, xy, X0, Y0, W0, bw);
 60                         else
 61                         #endif
 62                         for( ; x1 < bw; x1++ )
 63                         {
 64                             //計算W
 65                             double W = W0 + M[6]*x1;
 66                             W = W ? 1./W : 0;
 67                             //從齊次座標變到增廣矢量，注意如今是分塊操做，因此纔有+ M[0]*x1這一項
 68                             //把結果控制在(INT_MIN, INT_MAX)中，因爲下面會有把double轉爲int
 69                             double fX = std::max((double)INT_MIN, std::min((double)INT_MAX, (X0 + M[0]*x1)*W));
 70                             double fY = std::max((double)INT_MIN, std::min((double)INT_MAX, (Y0 + M[3]*x1)*W));
 71                             int X = saturate_cast<int>(fX);
 72                             int Y = saturate_cast<int>(fY);
 73 
 74                             xy[x1*2] = saturate_cast<short>(X);
 75                             xy[x1*2+1] = saturate_cast<short>(Y);
 76                         }
 77                     }
 78                     //其餘插值方法
 79                     else
 80                     {
 81                         short* alpha = A + y1*bw;
 82                         x1 = 0;
 83 
 84                         //原理上同
 85                         #if CV_TRY_SSE4_1
 86                         if (pwarp_impl_sse4)
 87                             pwarp_impl_sse4->process(M, xy, alpha, X0, Y0, W0, bw);
 88                         else
 89                         #endif
 90                         for( ;x1 < bw; x1++ )
 91                         {
 92                             double W = W0 + M[6]*x1;
 93                             //INTER_TAB_SIZE在個人opencv版本就是32，INTER_BITS是5
 94                             //下面的代碼是先把值擴大32倍，再左移5位獲得一個值
 95                             W = W ? INTER_TAB_SIZE/W : 0;
 96                             double fX = std::max((double)INT_MIN, std::min((double)INT_MAX, (X0 + M[0]*x1)*W));
 97                             double fY = std::max((double)INT_MIN, std::min((double)INT_MAX, (Y0 + M[3]*x1)*W));
 98                             int X = saturate_cast<int>(fX);
 99                             int Y = saturate_cast<int>(fY);
100 
101                             xy[x1*2] = saturate_cast<short>(X >> INTER_BITS);
102                             xy[x1*2+1] = saturate_cast<short>(Y >> INTER_BITS);
103                             
104                             //alpha是remap裏面所使用的remapBilinear所要用到的一個參數，這裏就不往下看了
105                             alpha[x1] = (short)((Y & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE +
106                                                 (X & (INTER_TAB_SIZE-1)));
107                         }
108                     }
109                 }
110 
111                 if( interpolation == INTER_NEAREST )
112                     //直接運行重映射，直接把源點映射到目標位置的點
113                     remap( src, dpart, _XY, Mat(), interpolation, borderType, borderValue );
114                 else
115                 {
116                     Mat _matA(bh, bw, CV_16U, A);
117                     remap( src, dpart, _XY, _matA, interpolation, borderType, borderValue );
118                 }
119             }
120         }
121     }
122 
123 private:
124     Mat src;
125     Mat dst;
126     const double* M;
127     int interpolation, borderType;
128     Scalar borderValue;
129 };

 

OpenCV的仿射變換源碼和投影變換的相似，只是仿射變換的矩陣M的元素個數是6個

 

參考

http://www.cnblogs.com/houkai/p/6660272.html

https://en.wikipedia.org/wiki/Bilinear_interpolation

http://www.cnblogs.com/tiandsp/archive/2012/12/16/2820916.html