OpenCV的代码提取：warpPerspective函数的实现

实现代码warpPerspective.hpp：

1. // fbc_cv是免费软件，并且使用与OpenCV相同的许可证  
2. #ifndef FBC_CV_WARP_PERSPECTIVE_HPP_  
3. #define FBC_CV_WARP_PERSPECTIVE_HPP_  
4.   
5. /* reference: include/opencv2/imgproc.hpp 
6.               modules/imgproc/src/imgwarp.cpp 
7. */  
8.   
9. #include <typeinfo>  
10. #include "core/mat.hpp"  
11. #include "core/invert.hpp"  
12. #include "imgproc.hpp"  
13. #include "remap.hpp"  
14.   
15. namespace fbc {  
16.   
17. // 从四对对应点计算透视变换  
18. FBC_EXPORTS int getPerspectiveTransform(const Point2f src1[], const Point2f src2[], Mat_<double, 1>& dst);  
19.   
20. // 对图像应用透视变换  
21. // 该功能无法正常工作  
22. // 支持类型：uchar / float  
23. /*  
24. \f[\texttt{ dst } (x, y) = \texttt{ src } \left(\frac{ M_{ 11 } x + M_{ 12 } y + M_{ 13 } }{M_{ 31 } x + M_{ 32 } y + M_{ 33 }}, 
25.     \frac{ M_{ 21 } x + M_{ 22 } y + M_{ 23 } }{M_{ 31 } x + M_{ 32 } y + M_{ 33 }} \right)\f] 
26. */  
27. template<typename _Tp1, typename _Tp2, int chs1, int chs2>  
28. int warpPerspective(const Mat_<_Tp1, chs1>& src, Mat_<_Tp1, chs1>& dst, const Mat_<_Tp2, chs2>& M_,  
29.     int flags = INTER_LINEAR, int borderMode = BORDER_CONSTANT, const Scalar& borderValue = Scalar())  
30. {  
31.     FBC_Assert(src.data != NULL && dst.data != NULL && M_.data != NULL);  
32.     FBC_Assert(src.cols > 0 && src.rows > 0 && dst.cols > 0 && dst.rows > 0);  
33.     FBC_Assert(src.data != dst.data);  
34.     FBC_Assert(typeid(double).name() == typeid(_Tp2).name() && M_.rows == 3 && M_.cols == 3);  
35.     FBC_Assert((typeid(uchar).name() == typeid(_Tp1).name()) || (typeid(float).name() == typeid(_Tp1).name())); // uchar/float  
36.   
37.     double M[9];  
38.     Mat_<double, 1> matM(3, 3, M);  
39.   
40.     int interpolation = flags & INTER_MAX;  
41.     if (interpolation == INTER_AREA)  
42.         interpolation = INTER_LINEAR;  
43.   
44.     if (!(flags & WARP_INVERSE_MAP))  
45.         invert(M_, matM);  
46.   
47.     Range range(0, dst.rows);  
48.   
49.     const int BLOCK_SZ = 32;  
50.     short XY[BLOCK_SZ*BLOCK_SZ * 2], A[BLOCK_SZ*BLOCK_SZ];  
51.     int x, y, x1, y1, width = dst.cols, height = dst.rows;  
52.   
53.     int bh0 = std::min(BLOCK_SZ / 2, height);  
54.     int bw0 = std::min(BLOCK_SZ*BLOCK_SZ / bh0, width);  
55.     bh0 = std::min(BLOCK_SZ*BLOCK_SZ / bw0, height);  
56.   
57.     for (y = range.start; y < range.end; y += bh0) {  
58.         for (x = 0; x < width; x += bw0) {  
59.             int bw = std::min(bw0, width - x);  
60.             int bh = std::min(bh0, range.end - y); // height  
61.   
62.             Mat_<short, 2> _XY(bh, bw, XY), matA;  
63.             Mat_<_Tp1, chs1> dpart;  
64.             dst.getROI(dpart, Rect(x, y, bw, bh));  
65.   
66.             for (y1 = 0; y1 < bh; y1++) {  
67.                 short* xy = XY + y1*bw * 2;  
68.                 double X0 = M[0] * x + M[1] * (y + y1) + M[2];  
69.                 double Y0 = M[3] * x + M[4] * (y + y1) + M[5];  
70.                 double W0 = M[6] * x + M[7] * (y + y1) + M[8];  
71.   
72.                 if (interpolation == INTER_NEAREST) {  
73.                     x1 = 0;  
74.                     for (; x1 < bw; x1++) {  
75.                         double W = W0 + M[6] * x1;  
76.                         W = W ? 1. / W : 0;  
77.                         double fX = std::max((double)INT_MIN, std::min((double)INT_MAX, (X0 + M[0] * x1)*W));  
78.                         double fY = std::max((double)INT_MIN, std::min((double)INT_MAX, (Y0 + M[3] * x1)*W));  
79.                         int X = saturate_cast<int>(fX);  
80.                         int Y = saturate_cast<int>(fY);  
81.   
82.                         xy[x1 * 2] = saturate_cast<short>(X);  
83.                         xy[x1 * 2 + 1] = saturate_cast<short>(Y);  
84.                     }  
85.                 } else {  
86.                     short* alpha = A + y1*bw;  
87.                     x1 = 0;  
88.                     for (; x1 < bw; x1++) {  
89.                         double W = W0 + M[6] * x1;  
90.                         W = W ? INTER_TAB_SIZE / W : 0;  
91.                         double fX = std::max((double)INT_MIN, std::min((double)INT_MAX, (X0 + M[0] * x1)*W));  
92.                         double fY = std::max((double)INT_MIN, std::min((double)INT_MAX, (Y0 + M[3] * x1)*W));  
93.                         int X = saturate_cast<int>(fX);  
94.                         int Y = saturate_cast<int>(fY);  
95.   
96.                         xy[x1 * 2] = saturate_cast<short>(X >> INTER_BITS);  
97.                         xy[x1 * 2 + 1] = saturate_cast<short>(Y >> INTER_BITS);  
98.                         alpha[x1] = (short)((Y & (INTER_TAB_SIZE - 1))*INTER_TAB_SIZE + (X & (INTER_TAB_SIZE - 1)));  
99.                     }  
100.                 }  
101.             }  
102.   
103.             if (interpolation == INTER_NEAREST) {  
104.                 remap(src, dpart, _XY, Mat_<float, 1>(), interpolation, borderMode, borderValue);  
105.             } else {  
106.                 Mat_<ushort, 1> _matA(bh, bw, A);  
107.                 remap(src, dpart, _XY, _matA, interpolation, borderMode, borderValue);  
108.             }  
109.         }  
110.     }  
111.   
112.     return 0;  
113. }  
114.   
115. } // 命名空间fbc  
116.   
117. #endif // FBC_CV_WARP_PERSPECTIVE_HPP_