图像放大并进行BiCubic插值 Matlab/C++代码


BiCubic 双三次插值


BiCubic插值原理:

双三次插值又称立方卷积插值。三次卷积插值是一种更加复杂的插值方式。该算法利用待采样点周围16个点的灰度值作三次插值,不仅考虑到4 个直接相邻点的灰度影响,而且考虑到各邻点间灰度值变化率的影响。三次运算可以得到更接近高分辨率图像的放大效果,但也导致了运算量的急剧增加。这种算法需要选取插值基函数来拟合数据,其最常用的插值基函数如图1所示,本次实验采用如图所示函数作为基函数。

构造BiCubic函数:

rbf插值算法python birkhoff插值_基函数

其中,a取-0.5.

BiCubic函数具有如下形状:

rbf插值算法python birkhoff插值_rbf插值算法python_02

[source:  R. Keys, (1981). "Cubic convolution interpolation for digital image processing". IEEE Transactions on Signal Processing, Acoustics, Speech, and Signal Processing 29 (6): 1153–1160.]

对待插值的像素点(x,y)(x和y可以为浮点数),取其附近的4x4邻域点(xi,yj), i,j = 0,1,2,3。按如下公式进行插值计算:

rbf插值算法python birkhoff插值_rbf插值算法python_03

Matlab实现代码:



[plain]  view plain copy


1. %双三次插值具体实现  
2. clc,clear;  
3. fff=imread('E:\Documents\BUPT\DIP\图片\lena.bmp');   
4. ff =rgb2gray(fff);%转化为灰度图像  
5. [mm,nn]=size(ff);               %将图像隔行隔列抽取元素,得到缩小的图像f  
6. m=mm/2;  
7. n=nn/2;  
8. f =zeros(m,n);  
9. for i=1:m  
10.    for j=1:n  
11.      f(i,j)=ff(2*i,2*j);  
12.    end  
13. end



[plain]  view plain copy



    1. k=5;                       %设置放大倍数  
    2. bijiao1 =imresize(f,k,'bilinear');%双线性插值结果比较  
    3. bijiao =uint8(bijiao1);



    [plain]  view plain copy


    1. a=f(1,:);  
    2. c=f(m,:);             %将待插值图像矩阵前后各扩展两行两列,共扩展四行四列  
    3. b=[f(1,1),f(1,1),f(:,1)',f(m,1),f(m,1)];  
    4. d=[f(1,n),f(1,n),f(:,n)',f(m,n),f(m,n)];  
    5. a1=[a;a;f;c;c];  
    6. b1=[b;b;a1';d;d];  
    7. ffff=b1';  
    8. f1=double(ffff);  
    9. g1 =zeros(k*m,k*n);  
    10. fori=1:k*m                 %利用双三次插值公式对新图象所有像素赋值  
    11.    u=rem(i,k)/k;  
    12. i1=floor(i/k)+2;  
    13.    A=[sw(1+u) sw(u) sw(1-u) sw(2-u)];    
    14.   for j=1:k*n  
    15.      v=rem(j,k)/k;  
    16. j1=floor(j/k)+2;  
    17.      C=[sw(1+v);sw(v);sw(1-v);sw(2-v)];  
    18.      B=[f1(i1-1,j1-1) f1(i1-1,j1) f1(i1-1,j1+1)f1(i1-1,j1+2)  
    19.        f1(i1,j1-1)   f1(i1,j1)  f1(i1,j1+1)   f1(i1,j1+2)  
    20.        f1(i1+1,j1-1)   f1(i1+1,j1) f1(i1+1,j1+1) f1(i1+1,j1+2)  
    21.        f1(i1+2,j1-1) f1(i1+2,j1) f1(i1+2,j1+1)f1(i1+2,j1+2)];  
    22.      g1(i,j)=(A*B*C);  
    23.    end  
    24. end  
    25. g=uint8(g1);


    [plain]  view plain copy


    1. imshow(uint8(f));title('缩小的图像');             %显示缩小的图像  
    2. figure,imshow(ff);title('原图');               %显示原图像  
    3. figure,imshow(g);title('双三次插值放大的图像');     %显示插值后的图像  
    4. figure,imshow(bijiao);title('双线性插值放大结果');     %显示插值后的图像   
    5. mse=0;  
    6. ff=double(ff);  
    7. g=double(g);              
    8. ff2=fftshift(fft2(ff));   %计算原图像和插值图像的傅立叶幅度谱                              
    9. g2=fftshift(fft2(g));  
    10. figure,subplot(1,2,1),imshow(log(abs(ff2)),[8,10]);title('原图像的傅立叶幅度谱');  
    11. subplot(1,2,2),imshow(log(abs(g2)),[8,10]);title('双三次插值图像的傅立叶幅度谱');

    [plain]  view plain copy



      1. 基函数代码:  
      2. functionA=sw(w1)  
      3. w=abs(w1);  
      4. ifw<1&&w>=0  
      5.    A=1-2*w^2+w^3;  
      6. elseifw>=1&&w<2  
      7.    A=4-8*w+5*w^2-w^3;  
      8. else  
      9.   A=0;  
      10. end




      C++实现代码:



      [html]  view plain copy


      1. <pre code_snippet_id="307700" snippet_file_name="blog_20140423_6_4966111" name="code" class="cpp" style="font-family: Arial;">#include "opencv2/imgproc/imgproc.hpp"  
      2. #include "opencv2/highgui/highgui.hpp"  
      3. #include <iostream>  
      4. #include <cmath>  
      5. #include <fstream>  
      6. using namespace cv;  
      7. using namespace std;  
      8. #define PI 3.14159265  
      9. float BiCubicPoly(float x);  
      10. void MyScaleBiCubicInter(Mat& src, Mat& dst, float TransMat[3][3]);  
      11. /**  
      12.  * @function main  
      13.  */  
      14. int main( int argc, char** argv )  
      15. {  
      16.   // load image  
      17. imageName = "images/Lenna_256.png";  
      18.   Mat image;  
      19. image = imread(imageName,1);</pre><pre code_snippet_id="307700" snippet_file_name="blog_20140423_7_964140" name="code" class="cpp" style="font-family: Arial;">  if(!image.data)  
      20.   {  
      21. << "No image data" << endl;  
      22.       return -1;  
      23.   }  
      24.   // show image  
      25.   namedWindow("image", CV_WINDOW_AUTOSIZE);  
      26.   imshow("image", image);    
      27.   Mat dst;  
      28. </pre><pre code_snippet_id="307700" snippet_file_name="blog_20140423_8_898884" name="code" class="cpp" style="font-family: Arial;">  MyScaleBiCubicInter(image, dst, transMat);  
      29.   namedWindow("out_image", CV_WINDOW_AUTOSIZE);  
      30.   imshow("out_image", dst);  
      31.   imwrite("Lenna_scale_biCubic2.jpg", dst);  
      32.   waitKey(0);  
      33.   return 0;  
      34. }  
      35. float BiCubicPoly(float x)  
      36. {  
      37. abs_x = abs(x);  
      38. a = -0.5;  
      39. <= 1.0 )  
      40.     {  
      41.         return (a+2)*pow(abs_x,3) - (a+3)*pow(abs_x,2) + 1;  
      42.     }  
      43. < 2.0 )  
      44.     {  
      45.         return a*pow(abs_x,3) - 5*a*pow(abs_x,2) + 8*a*abs_x - 4*a;  
      46.     }  
      47.     else  
      48.         return 0.0;  
      49. }  
      50.   
      51. void MyScaleBiCubicInter(Mat& src, Mat& dst, float TransMat[3][3])  
      52. {  
      53.     CV_Assert(src.data);  
      54.     CV_Assert(src.depth() != sizeof(uchar));  
      55.       
      56.     // calculate margin point of dst image  
      57. left =  0;  
      58. right =  0;  
      59. top =  0;  
      60. down =  0;  
      61.   
      62. x = src.cols * 1.0f;  
      63. y = 0.0f;  
      64. u1 = x * TransMat[0][0] + y * TransMat[0][1];  
      65. v1 = x * TransMat[1][0] + y * TransMat[1][1];  
      66. x = src.cols * 1.0f;  
      67. y = src.rows * 1.0f;  
      68. u2 = x * TransMat[0][0] + y * TransMat[0][1];  
      69. v2 = x * TransMat[1][0] + y * TransMat[1][1];  
      70. x = 0.0f;  
      71. y = src.rows * 1.0f;  
      72. u3 = x * TransMat[0][0] + y * TransMat[0][1];  
      73. v3 = x * TransMat[1][0] + y * TransMat[1][1];  
      74.   
      75. left =  min( min( min(0.0f,u1), u2 ), u3);  
      76. right =  max( max( max(0.0f,u1), u2 ), u3);  
      77. top =  min( min( min(0.0f,v1), v2 ), v3);  
      78. down =  max( max( max(0.0f,v1), v2 ), v3);  
      79.   
      80.     // create dst image  
      81.     dst.create(int(abs(right-left)), int(abs(down-top)), src.type());     
      82.   
      83.     CV_Assert( dst.channels() == src.channels() );  
      84. channels = dst.channels();  
      85.   
      86.     int i,j;  
      87.     uchar* p;  
      88.     uchar* q0;  
      89.     uchar* q1;  
      90.     uchar* q2;  
      91.     uchar* q3;  
      92. i = 0; i < dst.rows; ++i)  
      93.     {  
      94. p = dst.ptr<uchar>(i);  
      95. j = 0; j < dst.cols; ++j)  
      96.         {  
      97.             //   
      98. x = (j+left)/TransMat[0][0]  ;   
      99. y = (i+top)/TransMat[1][1] ;  
      100.   
      101. x0 = int(x) - 1;  
      102. y0 = int(y) - 1;  
      103. x1 = int(x);  
      104. y1 = int(y);  
      105. x2 = int(x) + 1;  
      106. y2 = int(y) + 1;  
      107. x3 = int(x) + 2;  
      108. y3 = int(y) + 2;  
      109.   
      110. >= 0) && (x3 < src.cols) && (y0 >= 0) && (y3 < src.rows) )   
      111.             {  
      112. q0 = src.ptr<uchar>(y0);  
      113. q1 = src.ptr<uchar>(y1);  
      114. q2 = src.ptr<uchar>(y2);  
      115. q3 = src.ptr<uchar>(y3);  
      116.                   
      117. dist_x0 = BiCubicPoly(x-x0);  
      118. dist_x1 = BiCubicPoly(x-x1);  
      119. dist_x2 = BiCubicPoly(x-x2);  
      120. dist_x3 = BiCubicPoly(x-x3);  
      121. dist_y0 = BiCubicPoly(y-y0);  
      122. dist_y1 = BiCubicPoly(y-y1);  
      123. dist_y2 = BiCubicPoly(y-y2);  
      124. dist_y3 = BiCubicPoly(y-y3);  
      125.   
      126. dist_x0y0 = dist_x0 * dist_y0;  
      127. dist_x0y1 = dist_x0 * dist_y1;  
      128. dist_x0y2 = dist_x0 * dist_y2;  
      129. dist_x0y3 = dist_x0 * dist_y3;  
      130. dist_x1y0 = dist_x1 * dist_y0;  
      131. dist_x1y1 = dist_x1 * dist_y1;  
      132. dist_x1y2 = dist_x1 * dist_y2;  
      133. dist_x1y3 = dist_x1 * dist_y3;  
      134. dist_x2y0 = dist_x2 * dist_y0;  
      135. dist_x2y1 = dist_x2 * dist_y1;  
      136. dist_x2y2 = dist_x2 * dist_y2;  
      137. dist_x2y3 = dist_x2 * dist_y3;  
      138. dist_x3y0 = dist_x3 * dist_y0;  
      139. dist_x3y1 = dist_x3 * dist_y1;  
      140. dist_x3y2 = dist_x3 * dist_y2;  
      141. dist_x3y3 = dist_x3 * dist_y3;  
      142.                   
      143.                 switch(channels)  
      144.                 {  
      145.                     case 1:  
      146.                         {  
      147.                             break;  
      148.                         }  
      149.                     case 3:  
      150.                         {  
      151.                             p[3*j] =    (uchar)(q0[3*x0] * dist_x0y0 +  
      152.                                                 q1[3*x0] * dist_x0y1 +  
      153.                                                 q2[3*x0] * dist_x0y2 +  
      154.                                                 q3[3*x0] * dist_x0y3 +  
      155.                                                 q0[3*x1] * dist_x1y0 +  
      156.                                                 q1[3*x1] * dist_x1y1 +  
      157.                                                 q2[3*x1] * dist_x1y2 +  
      158.                                                 q3[3*x1] * dist_x1y3 +  
      159.                                                 q0[3*x2] * dist_x2y0 +  
      160.                                                 q1[3*x2] * dist_x2y1 +  
      161.                                                 q2[3*x2] * dist_x2y2 +  
      162.                                                 q3[3*x2] * dist_x2y3 +  
      163.                                                 q0[3*x3] * dist_x3y0 +  
      164.                                                 q1[3*x3] * dist_x3y1 +  
      165.                                                 q2[3*x3] * dist_x3y2 +  
      166.                                                 q3[3*x3] * dist_x3y3 ) ;  
      167.   
      168.                             p[3*j+1] =  (uchar)(q0[3*x0+1] * dist_x0y0 +  
      169.                                                 q1[3*x0+1] * dist_x0y1 +  
      170.                                                 q2[3*x0+1] * dist_x0y2 +  
      171.                                                 q3[3*x0+1] * dist_x0y3 +  
      172.                                                 q0[3*x1+1] * dist_x1y0 +  
      173.                                                 q1[3*x1+1] * dist_x1y1 +  
      174.                                                 q2[3*x1+1] * dist_x1y2 +  
      175.                                                 q3[3*x1+1] * dist_x1y3 +  
      176.                                                 q0[3*x2+1] * dist_x2y0 +  
      177.                                                 q1[3*x2+1] * dist_x2y1 +  
      178.                                                 q2[3*x2+1] * dist_x2y2 +  
      179.                                                 q3[3*x2+1] * dist_x2y3 +  
      180.                                                 q0[3*x3+1] * dist_x3y0 +  
      181.                                                 q1[3*x3+1] * dist_x3y1 +  
      182.                                                 q2[3*x3+1] * dist_x3y2 +  
      183.                                                 q3[3*x3+1] * dist_x3y3 ) ;  
      184.   
      185.                             p[3*j+2] =  (uchar)(q0[3*x0+2] * dist_x0y0 +  
      186.                                                 q1[3*x0+2] * dist_x0y1 +  
      187.                                                 q2[3*x0+2] * dist_x0y2 +  
      188.                                                 q3[3*x0+2] * dist_x0y3 +  
      189.                                                 q0[3*x1+2] * dist_x1y0 +  
      190.                                                 q1[3*x1+2] * dist_x1y1 +  
      191.                                                 q2[3*x1+2] * dist_x1y2 +  
      192.                                                 q3[3*x1+2] * dist_x1y3 +  
      193.                                                 q0[3*x2+2] * dist_x2y0 +  
      194.                                                 q1[3*x2+2] * dist_x2y1 +  
      195.                                                 q2[3*x2+2] * dist_x2y2 +  
      196.                                                 q3[3*x2+2] * dist_x2y3 +  
      197.                                                 q0[3*x3+2] * dist_x3y0 +  
      198.                                                 q1[3*x3+2] * dist_x3y1 +  
      199.                                                 q2[3*x3+2] * dist_x3y2 +  
      200.                                                 q3[3*x3+2] * dist_x3y3 ) ;  
      201.   
      202. thre = 198.0f;  
      203. > thre) || (abs(p[3*j+1]-q1[3*x1+1]) > thre) ||  
      204. > thre) )  
      205.                             {  
      206.                                 p[3*j] = q1[3*x1];  
      207.                                 p[3*j+1] = q1[3*x1+1];  
      208.                                 p[3*j+2] = q1[3*x1+2];  
      209.                             }                     
      210.                             break;  
      211.                         }  
      212.                 }  
      213.             }  
      214.         }  
      215.     }  
      216. }</pre><pre code_snippet_id="307700" snippet_file_name="blog_20140423_9_8881779" name="code" class="cpp"><span style="font-family:Arial;">参考:</span><h2 style="margin: 0px 20px 0px 0px; padding: 0px; display: inline-block; line-height: 30px; background-color: rgb(204, 232, 207);"><a target="_blank" href="" style="text-decoration: none;"><span style="font-family:SimSun;font-size:18px;">lichengyu的专栏</span></a></h2></pre><p></p>  
      217. <pre></pre>  
      218. <p></p>  
      219.

      BiCubic插值原理:

      双三次插值又称立方卷积插值。三次卷积插值是一种更加复杂的插值方式。该算法利用待采样点周围16个点的灰度值作三次插值,不仅考虑到4 个直接相邻点的灰度影响,而且考虑到各邻点间灰度值变化率的影响。三次运算可以得到更接近高分辨率图像的放大效果,但也导致了运算量的急剧增加。这种算法需要选取插值基函数来拟合数据,其最常用的插值基函数如图1所示,本次实验采用如图所示函数作为基函数。

      构造BiCubic函数:

      rbf插值算法python birkhoff插值_基函数

      其中,a取-0.5.

      BiCubic函数具有如下形状:

      rbf插值算法python birkhoff插值_rbf插值算法python_02

      [source:  R. Keys, (1981). "Cubic convolution interpolation for digital image processing". IEEE Transactions on Signal Processing, Acoustics, Speech, and Signal Processing 29 (6): 1153–1160.]

      对待插值的像素点(x,y)(x和y可以为浮点数),取其附近的4x4邻域点(xi,yj), i,j = 0,1,2,3。按如下公式进行插值计算:

      rbf插值算法python birkhoff插值_rbf插值算法python_03

      Matlab实现代码:



      [plain]  view plain copy


      1. %双三次插值具体实现  
      2. clc,clear;  
      3. fff=imread('E:\Documents\BUPT\DIP\图片\lena.bmp');   
      4. ff =rgb2gray(fff);%转化为灰度图像  
      5. [mm,nn]=size(ff);               %将图像隔行隔列抽取元素,得到缩小的图像f  
      6. m=mm/2;  
      7. n=nn/2;  
      8. f =zeros(m,n);  
      9. for i=1:m  
      10.    for j=1:n  
      11.      f(i,j)=ff(2*i,2*j);  
      12.    end  
      13. end


      [plain]  view plain copy



        1. k=5;                       %设置放大倍数  
        2. bijiao1 =imresize(f,k,'bilinear');%双线性插值结果比较  
        3. bijiao =uint8(bijiao1);


        [plain]  view plain copy


        1. a=f(1,:);  
        2. c=f(m,:);             %将待插值图像矩阵前后各扩展两行两列,共扩展四行四列  
        3. b=[f(1,1),f(1,1),f(:,1)',f(m,1),f(m,1)];  
        4. d=[f(1,n),f(1,n),f(:,n)',f(m,n),f(m,n)];  
        5. a1=[a;a;f;c;c];  
        6. b1=[b;b;a1';d;d];  
        7. ffff=b1';  
        8. f1=double(ffff);  
        9. g1 =zeros(k*m,k*n);  
        10. fori=1:k*m                 %利用双三次插值公式对新图象所有像素赋值  
        11.    u=rem(i,k)/k;  
        12. i1=floor(i/k)+2;  
        13.    A=[sw(1+u) sw(u) sw(1-u) sw(2-u)];    
        14.   for j=1:k*n  
        15.      v=rem(j,k)/k;  
        16. j1=floor(j/k)+2;  
        17.      C=[sw(1+v);sw(v);sw(1-v);sw(2-v)];  
        18.      B=[f1(i1-1,j1-1) f1(i1-1,j1) f1(i1-1,j1+1)f1(i1-1,j1+2)  
        19.        f1(i1,j1-1)   f1(i1,j1)  f1(i1,j1+1)   f1(i1,j1+2)  
        20.        f1(i1+1,j1-1)   f1(i1+1,j1) f1(i1+1,j1+1) f1(i1+1,j1+2)  
        21.        f1(i1+2,j1-1) f1(i1+2,j1) f1(i1+2,j1+1)f1(i1+2,j1+2)];  
        22.      g1(i,j)=(A*B*C);  
        23.    end  
        24. end  
        25. g=uint8(g1);



        [plain]  view plain copy



          1. imshow(uint8(f));title('缩小的图像');             %显示缩小的图像  
          2. figure,imshow(ff);title('原图');               %显示原图像  
          3. figure,imshow(g);title('双三次插值放大的图像');     %显示插值后的图像  
          4. figure,imshow(bijiao);title('双线性插值放大结果');     %显示插值后的图像   
          5. mse=0;  
          6. ff=double(ff);  
          7. g=double(g);              
          8. ff2=fftshift(fft2(ff));   %计算原图像和插值图像的傅立叶幅度谱                              
          9. g2=fftshift(fft2(g));  
          10. figure,subplot(1,2,1),imshow(log(abs(ff2)),[8,10]);title('原图像的傅立叶幅度谱');  
          11. subplot(1,2,2),imshow(log(abs(g2)),[8,10]);title('双三次插值图像的傅立叶幅度谱');



          [plain]  view plain copy



          1. 基函数代码:  
          1. functionA=sw(w1)  
          2. w=abs(w1);  
          3. ifw<1&&w>=0  
          4.    A=1-2*w^2+w^3;  
          5. elseifw>=1&&w<2  
          6.    A=4-8*w+5*w^2-w^3;  
          7. else  
          8.   A=0;  
          9. end




          C++实现代码:



          [html]  view plain copy



          1. <pre code_snippet_id="307700" snippet_file_name="blog_20140423_6_4966111" name="code" class="cpp" style="font-family: Arial;">#include "opencv2/imgproc/imgproc.hpp"  
          2. #include "opencv2/highgui/highgui.hpp"  
          3. #include <iostream>  
          4. #include <cmath>  
          5. #include <fstream>  
          6. using namespace cv;  
          7. using namespace std;  
          8. #define PI 3.14159265  
          9. float BiCubicPoly(float x);  
          10. void MyScaleBiCubicInter(Mat& src, Mat& dst, float TransMat[3][3]);  
          11. /**  
          12.  * @function main  
          13.  */  
          14. int main( int argc, char** argv )  
          15. {  
          16.   // load image  
          17. imageName = "images/Lenna_256.png";  
          18.   Mat image;  
          19. image = imread(imageName,1);</pre><pre code_snippet_id="307700" snippet_file_name="blog_20140423_7_964140" name="code" class="cpp" style="font-family: Arial;">  if(!image.data)  
          20.   {  
          21. << "No image data" << endl;  
          22.       return -1;  
          23.   }  
          24.   // show image  
          25.   namedWindow("image", CV_WINDOW_AUTOSIZE);  
          26.   imshow("image", image);    
          27.   Mat dst;  
          28. </pre><pre code_snippet_id="307700" snippet_file_name="blog_20140423_8_898884" name="code" class="cpp" style="font-family: Arial;">  MyScaleBiCubicInter(image, dst, transMat);  
          29.   namedWindow("out_image", CV_WINDOW_AUTOSIZE);  
          30.   imshow("out_image", dst);  
          31.   imwrite("Lenna_scale_biCubic2.jpg", dst);  
          32.   waitKey(0);  
          33.   return 0;  
          34. }  
          35. float BiCubicPoly(float x)  
          36. {  
          37. abs_x = abs(x);  
          38. a = -0.5;  
          39. <= 1.0 )  
          40.     {  
          41.         return (a+2)*pow(abs_x,3) - (a+3)*pow(abs_x,2) + 1;  
          42.     }  
          43. < 2.0 )  
          44.     {  
          45.         return a*pow(abs_x,3) - 5*a*pow(abs_x,2) + 8*a*abs_x - 4*a;  
          46.     }  
          47.     else  
          48.         return 0.0;  
          49. }  
          50.   
          51. void MyScaleBiCubicInter(Mat& src, Mat& dst, float TransMat[3][3])  
          52. {  
          53.     CV_Assert(src.data);  
          54.     CV_Assert(src.depth() != sizeof(uchar));  
          55.       
          56.     // calculate margin point of dst image  
          57. left =  0;  
          58. right =  0;  
          59. top =  0;  
          60. down =  0;  
          61.   
          62. x = src.cols * 1.0f;  
          63. y = 0.0f;  
          64. u1 = x * TransMat[0][0] + y * TransMat[0][1];  
          65. v1 = x * TransMat[1][0] + y * TransMat[1][1];  
          66. x = src.cols * 1.0f;  
          67. y = src.rows * 1.0f;  
          68. u2 = x * TransMat[0][0] + y * TransMat[0][1];  
          69. v2 = x * TransMat[1][0] + y * TransMat[1][1];  
          70. x = 0.0f;  
          71. y = src.rows * 1.0f;  
          72. u3 = x * TransMat[0][0] + y * TransMat[0][1];  
          73. v3 = x * TransMat[1][0] + y * TransMat[1][1];  
          74.   
          75. left =  min( min( min(0.0f,u1), u2 ), u3);  
          76. right =  max( max( max(0.0f,u1), u2 ), u3);  
          77. top =  min( min( min(0.0f,v1), v2 ), v3);  
          78. down =  max( max( max(0.0f,v1), v2 ), v3);  
          79.   
          80.     // create dst image  
          81.     dst.create(int(abs(right-left)), int(abs(down-top)), src.type());     
          82.   
          83.     CV_Assert( dst.channels() == src.channels() );  
          84. channels = dst.channels();  
          85.   
          86.     int i,j;  
          87.     uchar* p;  
          88.     uchar* q0;  
          89.     uchar* q1;  
          90.     uchar* q2;  
          91.     uchar* q3;  
          92. i = 0; i < dst.rows; ++i)  
          93.     {  
          94. p = dst.ptr<uchar>(i);  
          95. j = 0; j < dst.cols; ++j)  
          96.         {  
          97.             //   
          98. x = (j+left)/TransMat[0][0]  ;   
          99. y = (i+top)/TransMat[1][1] ;  
          100.   
          101. x0 = int(x) - 1;  
          102. y0 = int(y) - 1;  
          103. x1 = int(x);  
          104. y1 = int(y);  
          105. x2 = int(x) + 1;  
          106. y2 = int(y) + 1;  
          107. x3 = int(x) + 2;  
          108. y3 = int(y) + 2;  
          109.   
          110. >= 0) && (x3 < src.cols) && (y0 >= 0) && (y3 < src.rows) )   
          111.             {  
          112. q0 = src.ptr<uchar>(y0);  
          113. q1 = src.ptr<uchar>(y1);  
          114. q2 = src.ptr<uchar>(y2);  
          115. q3 = src.ptr<uchar>(y3);  
          116.                   
          117. dist_x0 = BiCubicPoly(x-x0);  
          118. dist_x1 = BiCubicPoly(x-x1);  
          119. dist_x2 = BiCubicPoly(x-x2);  
          120. dist_x3 = BiCubicPoly(x-x3);  
          121. dist_y0 = BiCubicPoly(y-y0);  
          122. dist_y1 = BiCubicPoly(y-y1);  
          123. dist_y2 = BiCubicPoly(y-y2);  
          124. dist_y3 = BiCubicPoly(y-y3);  
          125.   
          126. dist_x0y0 = dist_x0 * dist_y0;  
          127. dist_x0y1 = dist_x0 * dist_y1;  
          128. dist_x0y2 = dist_x0 * dist_y2;  
          129. dist_x0y3 = dist_x0 * dist_y3;  
          130. dist_x1y0 = dist_x1 * dist_y0;  
          131. dist_x1y1 = dist_x1 * dist_y1;  
          132. dist_x1y2 = dist_x1 * dist_y2;  
          133. dist_x1y3 = dist_x1 * dist_y3;  
          134. dist_x2y0 = dist_x2 * dist_y0;  
          135. dist_x2y1 = dist_x2 * dist_y1;  
          136. dist_x2y2 = dist_x2 * dist_y2;  
          137. dist_x2y3 = dist_x2 * dist_y3;  
          138. dist_x3y0 = dist_x3 * dist_y0;  
          139. dist_x3y1 = dist_x3 * dist_y1;  
          140. dist_x3y2 = dist_x3 * dist_y2;  
          141. dist_x3y3 = dist_x3 * dist_y3;  
          142.                   
          143.                 switch(channels)  
          144.                 {  
          145.                     case 1:  
          146.                         {  
          147.                             break;  
          148.                         }  
          149.                     case 3:  
          150.                         {  
          151.                             p[3*j] =    (uchar)(q0[3*x0] * dist_x0y0 +  
          152.                                                 q1[3*x0] * dist_x0y1 +  
          153.                                                 q2[3*x0] * dist_x0y2 +  
          154.                                                 q3[3*x0] * dist_x0y3 +  
          155.                                                 q0[3*x1] * dist_x1y0 +  
          156.                                                 q1[3*x1] * dist_x1y1 +  
          157.                                                 q2[3*x1] * dist_x1y2 +  
          158.                                                 q3[3*x1] * dist_x1y3 +  
          159.                                                 q0[3*x2] * dist_x2y0 +  
          160.                                                 q1[3*x2] * dist_x2y1 +  
          161.                                                 q2[3*x2] * dist_x2y2 +  
          162.                                                 q3[3*x2] * dist_x2y3 +  
          163.                                                 q0[3*x3] * dist_x3y0 +  
          164.                                                 q1[3*x3] * dist_x3y1 +  
          165.                                                 q2[3*x3] * dist_x3y2 +  
          166.                                                 q3[3*x3] * dist_x3y3 ) ;  
          167.   
          168.                             p[3*j+1] =  (uchar)(q0[3*x0+1] * dist_x0y0 +  
          169.                                                 q1[3*x0+1] * dist_x0y1 +  
          170.                                                 q2[3*x0+1] * dist_x0y2 +  
          171.                                                 q3[3*x0+1] * dist_x0y3 +  
          172.                                                 q0[3*x1+1] * dist_x1y0 +  
          173.                                                 q1[3*x1+1] * dist_x1y1 +  
          174.                                                 q2[3*x1+1] * dist_x1y2 +  
          175.                                                 q3[3*x1+1] * dist_x1y3 +  
          176.                                                 q0[3*x2+1] * dist_x2y0 +  
          177.                                                 q1[3*x2+1] * dist_x2y1 +  
          178.                                                 q2[3*x2+1] * dist_x2y2 +  
          179.                                                 q3[3*x2+1] * dist_x2y3 +  
          180.                                                 q0[3*x3+1] * dist_x3y0 +  
          181.                                                 q1[3*x3+1] * dist_x3y1 +  
          182.                                                 q2[3*x3+1] * dist_x3y2 +  
          183.                                                 q3[3*x3+1] * dist_x3y3 ) ;  
          184.   
          185.                             p[3*j+2] =  (uchar)(q0[3*x0+2] * dist_x0y0 +  
          186.                                                 q1[3*x0+2] * dist_x0y1 +  
          187.                                                 q2[3*x0+2] * dist_x0y2 +  
          188.                                                 q3[3*x0+2] * dist_x0y3 +  
          189.                                                 q0[3*x1+2] * dist_x1y0 +  
          190.                                                 q1[3*x1+2] * dist_x1y1 +  
          191.                                                 q2[3*x1+2] * dist_x1y2 +  
          192.                                                 q3[3*x1+2] * dist_x1y3 +  
          193.                                                 q0[3*x2+2] * dist_x2y0 +  
          194.                                                 q1[3*x2+2] * dist_x2y1 +  
          195.                                                 q2[3*x2+2] * dist_x2y2 +  
          196.                                                 q3[3*x2+2] * dist_x2y3 +  
          197.                                                 q0[3*x3+2] * dist_x3y0 +  
          198.                                                 q1[3*x3+2] * dist_x3y1 +  
          199.                                                 q2[3*x3+2] * dist_x3y2 +  
          200.                                                 q3[3*x3+2] * dist_x3y3 ) ;  
          201.   
          202. thre = 198.0f;  
          203. > thre) || (abs(p[3*j+1]-q1[3*x1+1]) > thre) ||  
          204. > thre) )  
          205.                             {  
          206.                                 p[3*j] = q1[3*x1];  
          207.                                 p[3*j+1] = q1[3*x1+1];  
          208.                                 p[3*j+2] = q1[3*x1+2];  
          209.                             }                     
          210.                             break;  
          211.                         }  
          212.                 }  
          213.             }  
          214.         }  
          215.     }  
          216. }</pre><pre code_snippet_id="307700" snippet_file_name="blog_20140423_9_8881779" name="code" class="cpp"><span style="font-family:Arial;">参考:</span><h2 style="margin: 0px 20px 0px 0px; padding: 0px; display: inline-block; line-height: 30px; background-color: rgb(204, 232, 207);"><a target="_blank" href="" style="text-decoration: none;"><span style="font-family:SimSun;font-size:18px;">lichengyu的专栏</span></a></h2></pre><p></p>  
          217. <pre></pre>  
          218. <p></p>  
          219.