Image Resize 双线性和双立方插值法

I.图像处理-双线性插值（转载）

图像的缩放很好理解,就是图像的放大和缩小。传统的绘画工具中,有一种叫做“放大尺”的绘画工具，画家常用它来放大图画。当然，在计算机上，我们不再需要用放大尺去放大或缩小图像了，把这个工作交给程序来完成就可以了。下面就来讲讲计算机怎么来放大缩小图象；在本文中，我们所说的图像都是指点阵图，也就是用一个像素矩阵来描述图像的方法，对于另一种图像：用函数来描述图像的矢量图，不在本文讨论之列。
越是简单的模型越适合用来举例子，我们就举个简单的图像：3X3 的256级灰度图，也就是高为3个象素，宽也是3个象素的图像，每个象素的取值可以是 0－255，代表该像素的亮度，255代表最亮，也就是白色，0代表最暗，即黑色。假如图像的象素矩阵如下图所示（这个原始图把它叫做源图，Source）：
234   38    22
67     44    12
89     65    63

这个矩阵中，元素坐标(x,y)是这样确定的，x从左到右，从0开始，y从上到下，也是从零开始，这是图象处理中最常用的坐标系，就是这样一个坐标：

  ----------------------＞X
  |
  |
  |
  |
  |
∨Y

如果想把这副图放大为 4X4大小的图像，那么该怎么做呢？那么第一步肯定想到的是先把4X4的矩阵先画出来再说，好了矩阵画出来了，如下所示，当然，矩阵的每个像素都是未知数，等待着我们去填充（这个将要被填充的图的叫做目标图,Destination）：
?        ?        ?       ?
?        ?        ?       ?
?        ?        ?       ?
?        ?        ?       ? 
               
        然后要往这个空的矩阵里面填值是从源图中来，先填写目标图最左上角的象素，坐标为（0，0），那么该坐标对应源图中的坐标可以由如下公式得出：                                      
srcX=dstX* (srcWidth/dstWidth) , srcY = dstY * (srcHeight/dstHeight)
套用公式，就可以找到对应的原图的坐标(0*(3/4),0*(3/4))=>(0*0.75,0*0.75)=>(0,0)
,找到了源图的对应坐标,就可以把源图中坐标为(0,0)处的234象素值填进去目标图的(0,0)这个位置了。

接下来,如法炮制,寻找目标图中坐标为(1,0)的象素对应源图中的坐标,套用公式:
(1*0.75,0*0.75)=>(0.75,0)
结果发现,得到的坐标里面竟然有小数,这可怎么办?计算机里的图像可是数字图像,象素就是最小单位了,象素的坐标都是整数,从来没有小数坐标。这时候采用的一种策略就是采用四舍五入的方法（也可以采用直接舍掉小数位的方法），把非整数坐标转换成整数，好，那么按照四舍五入的方法就得到坐标（1，0），完整的运算过程就是这样的：
(1*0.75,0*0.75)=>(0.75,0)=>(1,0)
那么就可以再填一个象素到目标矩阵中了，同样是把源图中坐标为(1,0)处的像素值38填入目标图中的坐标。
         
依次填完每个象素，一幅放大后的图像就诞生了，像素矩阵如下所示：
234    38     22     22  
67      44     12     12  
89      65     63     63  
89      65     63     63  

这种放大图像的方法叫做最临近插值算法，这是一种最基本、最简单的图像缩放算法，效果也是最不好的，放大后的图像有很严重的马赛克，缩小后的图像有很严重的失真；效果不好的根源就是其简单的最临近插值方法引入了严重的图像失真，比如，当由目标图的坐标反推得到的源图的的坐标是一个浮点数的时候，采用了四舍五入的方法，直接采用了和这个浮点数最接近的象素的值，这种方法是很不科学的，当推得坐标值为 0.75的时候，不应该就简单的取为1，既然是0.75，比1要小0.25 ，比0要大0.75 ,那么目标象素值其实应该根据这个源图中虚拟的点四周的四个真实的点来按照一定的规律计算出来的，这样才能达到更好的缩放效果。双线型内插值算法就是一种比较好的图像缩放算法，它充分的利用了源图中虚拟点四周的四个真实存在的像素值来共同决定目标图中的一个像素值，因此缩放效果比简单的最邻近插值要好很多。

双线性内插值算法描述如下:
　　对于一个目的像素，设置坐标通过反向变换得到的浮点坐标为(i+u,j+v) (其中i、j均为浮点坐标的整数部分，u、v为浮点坐标的小数部分，是取值[0,1)区间的浮点数)，则这个像素得值 f(i+u,j+v) 可由原图像中坐标为 (i,j)、(i+1,j)、(i,j+1)、(i+1,j+1)所对应的周围四个像素的值决定，即：

f(i+u,j+v) = (1-u)(1-v)f(i,j) + (1-u)vf(i,j+1) + u(1-v)f(i+1,j)+ uvf(i+1,j+1)                          公式1

其中f(i,j)表示源图像(i,j)处的的像素值，以此类推。

比如，象刚才的例子，现在假如目标图的象素坐标为（1，1），那么反推得到的对应于源图的坐标是（0.75 , 0.75）, 这其实只是一个概念上的虚拟象素,实际在源图中并不存在这样一个象素,那么目标图的象素（1，1）的取值不能够由这个虚拟象素来决定，而只能由源图的这四个象素共同决定：（0，0）（0，1）（1，0）（1，1），而由于（0.75,0.75）离（1，1）要更近一些，那么（1,1）所起的决定作用更大一些，这从公式1中的系数uv=0.75×0.75就可以体现出来，而（0.75,0.75）离（0，0）最远，所以（0，0）所起的决定作用就要小一些，公式中系数为(1-u)(1-v)=0.25×0.25也体现出了这一特点；

最邻近插值和双向性内插值缩放图片的效果对比：

    

 原始图片

     

最邻近插值放大图片

        

双线型内插值放大图片

 

常用的插值方法有：最近邻插值、双线性插值、三次卷积法。

 

在做数字图像处理时，经常会碰到小数象素坐标的取值问题，这时就需要依据邻近象素的值来对该坐标进行插值。比如：做地图投影转换，对目标图像的一个象素进行坐标变换到源图像上对应的点时，变换出来的对应的坐标是一个小数，再比如做图像的几何校正，也会碰到同样的问题。以下是对常用的三种数字图像插值方法进行介绍。

1、最邻近元法

　　这是最简单的一种插值方法，不需要计算，在待求象素的四邻象素中，将距离待求象素最近的邻象素灰度赋给待求象素。设i+u, j+v(i, j为正整数， u, v为大于零小于1的小数，下同)为待求象素坐标，则待求象素灰度的值 f(i+u, j+v)　如下图所示：

 

 

如果(i+u, j+v)落在A区，即u<0.5, v<0.5，则将左上角象素的灰度值赋给待求象素，同理，落在B区则赋予右上角的象素灰度值，落在C区则赋予左下角象素的灰度值，落在D区则赋予右下角象素的灰度值。

最邻近元法计算量较小，但可能会造成插值生成的图像灰度上的不连续，在灰度变化的地方可能出现明显的锯齿状。

 

2、双线性内插法

双线性内插法是利用待求象素四个邻象素的灰度在两个方向上作线性内插，如下图所示：

 

对于 (i, j+v)，f(i, j) 到 f(i, j+1) 的灰度变化为线性关系，则有：

f(i, j+v) = [f(i, j+1)- f(i, j)] * v + f(i, j)

同理对于 (i+1, j+v) 则有：

                  f(i+1,j+v) = [f(i+1, j+1) - f(i+1, j)] * v + f(i+1, j)

从f(i, j+v) 到 f(i+1, j+v) 的灰度变化也为线性关系，由此可推导出待求象素灰度的计算式如下：

                  f(i+u,j+v) = (1-u) * (1-v) * f(i, j) + (1-u) * v * f(i, j+1) + u * (1-v) * f(i+1, j)+ u * v * f(i+1, j+1)

双线性内插法的计算比最邻近点法复杂，计算量较大，但没有灰度不连续的缺点，结果基本令人满意。它具有低通滤波性质，使高频分量受损，图像轮廓可能会有一点模糊。

3、三次内插法

该方法利用三次多项式S(x)求逼近理论上最佳插值函数sin(x)/x, 其数学表达式为：

待求像素(x, y)的灰度值由其周围16个灰度值加权内插得到，如下图：

 

 

待求像素的灰度计算式如下：

 

f(x, y) = f(i+u, j+v) = ABC

 

其中: 

 

三次曲线插值方法计算量较大，但插值后的图像效果最好。

II.图像处理-双立方插值（转载）

双立方插值计算涉及到16个像素点，其中(i’, j’)表示待计算像素点在源图像中的包含

小数部分的像素坐标，dx表示X方向的小数坐标，dy表示Y方向的小数坐标。具体

可以看下图：

根据上述图示与双立方插值的数学表达式可以看出，双立方插值本质上图像16个像素点

权重卷积之和作为新的像素值。

其中R(x)表示插值表达式，可以根据需要选择的表达式不同。常见有基于三角取值、Bell

分布表达、B样条曲线表达式。

1. 基于三角形采样数学公式为

最简单的线性分布，代码实现如下：

[java] view plain copy

1. private double triangleInterpolation( double f )  
2. {  
3.     f = f / 2.0;  
4.     if( f < 0.0 )  
5.     {  
6.         return ( f + 1.0 );  
7.     }  
8.     else  
9.     {  
10.         return ( 1.0 - f );  
11.     }  
12. }  

2.基于Bell分布采样的数学公式如下：

Bell分布采样数学公式基于三次卷积计算实现。代码实现如下：

[java] view plain copy

1. private double bellInterpolation( double x )  
2. {  
3.     double f = ( x / 2.0 ) * 1.5;  
4.     if( f > -1.5 && f < -0.5 )  
5.     {  
6.         return( 0.5 * Math.pow(f + 1.5, 2.0));  
7.     }  
8.     else if( f > -0.5 && f < 0.5 )  
9.     {  
10.         return 3.0 / 4.0 - ( f * f );  
11.     }  
12.     else if( ( f > 0.5 && f < 1.5 ) )  
13.     {  
14.         return( 0.5 * Math.pow(f - 1.5, 2.0));  
15.     }  
16.     return 0.0;  
17. }  

3.基于B样条曲线采样的数学公式如下：

是一种基于多项式的四次卷积的采样计算，代码如下：

[java] view plain copy

1. private double bspLineInterpolation( double f )  
2. {  
3.     if( f < 0.0 )  
4.     {  
5.         f = -f;  
6.     }  
7.   
8.     if( f >= 0.0 && f <= 1.0 )  
9.     {  
10.         return ( 2.0 / 3.0 ) + ( 0.5 ) * ( f* f * f ) - (f*f);  
11.     }  
12.     else if( f > 1.0 && f <= 2.0 )  
13.     {  
14.         return 1.0 / 6.0 * Math.pow( ( 2.0 - f  ), 3.0 );  
15.     }  
16.     return 1.0;  
17. }  

实现图像双立方插值的完整源代码如下：

[java] view plain copy

1. package com.gloomyfish.zoom.study;  
2.   
3. import java.awt.image.BufferedImage;  
4. import java.awt.image.ColorModel;  
5.   
6. import com.gloomyfish.filter.study.AbstractBufferedImageOp;  
7.   
8. public class BicubicInterpolationFilter extends AbstractBufferedImageOp  {  
9.     public final static int TRIANGLE__INTERPOLATION = 1;  
10.     public final static int BELL__INTERPOLATION = 2;  
11.     public final static int BSPLINE__INTERPOLATION = 4;  
12.     public final static int CATMULLROOM__INTERPOLATION = 8;  
13.     public final static double B = 0.0;  
14.     public final static double C = 0.5; // constant  
15.     private int destH; // zoom height  
16.     private int destW; // zoom width  
17.     private int type;  
18.     public BicubicInterpolationFilter()  
19.     {  
20.         this.type = BSPLINE__INTERPOLATION;  
21.     }  
22.     public void setType(int type) {  
23.         this.type = type;  
24.     }  
25.     public void setDestHeight(int destH) {  
26.         this.destH = destH;  
27.     }  
28.   
29.     public void setDestWidth(int destW) {  
30.         this.destW = destW;  
31.     }  
32.       
33.     private double bellInterpolation( double x )  
34.     {  
35.         double f = ( x / 2.0 ) * 1.5;  
36.         if( f > -1.5 && f < -0.5 )  
37.         {  
38.             return( 0.5 * Math.pow(f + 1.5, 2.0));  
39.         }  
40.         else if( f > -0.5 && f < 0.5 )  
41.         {  
42.             return 3.0 / 4.0 - ( f * f );  
43.         }  
44.         else if( ( f > 0.5 && f < 1.5 ) )  
45.         {  
46.             return( 0.5 * Math.pow(f - 1.5, 2.0));  
47.         }  
48.         return 0.0;  
49.     }  
50.       
51.     private double bspLineInterpolation( double f )  
52.     {  
53.         if( f < 0.0 )  
54.         {  
55.             f = -f;  
56.         }  
57.   
58.         if( f >= 0.0 && f <= 1.0 )  
59.         {  
60.             return ( 2.0 / 3.0 ) + ( 0.5 ) * ( f* f * f ) - (f*f);  
61.         }  
62.         else if( f > 1.0 && f <= 2.0 )  
63.         {  
64.             return 1.0 / 6.0 * Math.pow( ( 2.0 - f  ), 3.0 );  
65.         }  
66.         return 1.0;  
67.     }  
68.       
69.     private double triangleInterpolation( double f )  
70.     {  
71.         f = f / 2.0;  
72.         if( f < 0.0 )  
73.         {  
74.             return ( f + 1.0 );  
75.         }  
76.         else  
77.         {  
78.             return ( 1.0 - f );  
79.         }  
80.     }  
81.       
82.     private double CatMullRomInterpolation( double f )  
83.     {  
84.         if( f < 0.0 )  
85.         {  
86.             f = Math.abs(f);  
87.         }  
88.         if( f < 1.0 )  
89.         {  
90.             return ( ( 12 - 9 * B - 6 * C ) * ( f * f * f ) +  
91.                 ( -18 + 12 * B + 6 *C ) * ( f * f ) +  
92.                 ( 6 - 2 * B ) ) / 6.0;  
93.         }  
94.         else if( f >= 1.0 && f < 2.0 )  
95.         {  
96.             return ( ( -B - 6 * C ) * ( f * f * f )  
97.                 + ( 6 * B + 30 * C ) * ( f *f ) +  
98.                 ( - ( 12 * B ) - 48 * C  ) * f +  
99.                 8 * B + 24 * C)/ 6.0;  
100.         }  
101.         else  
102.         {  
103.             return 0.0;  
104.         }  
105.     }   
106.   
107.     @Override  
108.     public BufferedImage filter(BufferedImage src, BufferedImage dest) {  
109.         int width = src.getWidth();  
110.         int height = src.getHeight();  
111.   
112.         if (dest == null)  
113.             dest = createCompatibleDestImage(src, null);  
114.   
115.         int[] inPixels = new int[width * height];  
116.         int[] outPixels = new int[destH * destW];  
117.         getRGB(src, 0, 0, width, height, inPixels);  
118.         float rowRatio = ((float) height) / ((float) destH);  
119.         float colRatio = ((float) width) / ((float) destW);  
120.         int index = 0;  
121.         for (int row = 0; row < destH; row++) {  
122.             int ta = 0, tr = 0, tg = 0, tb = 0;  
123.             double srcRow = ((float) row) * rowRatio;  
124.             // 获取整数部分坐标 row Index  
125.             double j = Math.floor(srcRow);  
126.             // 获取行的小数部分坐标  
127.             double t = srcRow - j;  
128.             for (int col = 0; col < destW; col++) {  
129.                 double srcCol = ((float) col) * colRatio;  
130.                 // 获取整数部分坐标 column Index  
131.                 double k = Math.floor(srcCol);  
132.                 // 获取列的小数部分坐标  
133.                 double u = srcCol - k;  
134.                 double[] rgbData = new double[3];  
135.                 double rgbCoffeData = 0.0;  
136.                 for(int m=-1; m<3; m++)  
137.                 {  
138.                     for(int n=-1; n<3; n++)  
139.                     {  
140.                         int[] rgb = getPixel(j+m, k+n, width, height, inPixels);  
141.                         double f1 = 0.0d;  
142.                         double f2 = 0.0d;  
143.                         if(type == TRIANGLE__INTERPOLATION)  
144.                         {  
145.                             f1  = triangleInterpolation( ((double) m ) - t );  
146.                             f2 = triangleInterpolation ( -(( (double) n ) - u ) );    
147.                         }  
148.                         else if(type == BELL__INTERPOLATION)  
149.                         {  
150.                             f1  = bellInterpolation( ((double) m ) - t );  
151.                             f2 = bellInterpolation ( -(( (double) n ) - u ) );    
152.                         }  
153.                         else if(type == BSPLINE__INTERPOLATION)  
154.                         {  
155.                             f1  = bspLineInterpolation( ((double) m ) - t );  
156.                             f2 = bspLineInterpolation ( -(( (double) n ) - u ) );     
157.                         }  
158.                         else  
159.                         {  
160.                             f1  = CatMullRomInterpolation( ((double) m ) - t );  
161.                             f2 = CatMullRomInterpolation ( -(( (double) n ) - u ) );                              
162.                         }  
163.                         // sum of weight  
164.                         rgbCoffeData += f2*f1;  
165.                         // sum of the RGB values  
166.                         rgbData[0] += rgb[0] * f2 * f1;  
167.                         rgbData[1] += rgb[1] * f2 * f1;  
168.                         rgbData[2] += rgb[2] * f2 * f1;  
169.                     }  
170.                 }  
171.                 ta = 255;  
172.                 // get Red/green/blue value for sample pixel  
173.                 tr = (int) (rgbData[0]/rgbCoffeData);  
174.                 tg = (int) (rgbData[1]/rgbCoffeData);  
175.                 tb = (int) (rgbData[2]/rgbCoffeData);  
176.                 index = row * destW + col;  
177.                 outPixels[index] = (ta << 24) | (clamp(tr) << 16)  
178.                         | (clamp(tg) << 8) | clamp(tb);  
179.             }  
180.         }  
181.         setRGB(dest, 0, 0, destW, destH, outPixels);  
182.         return dest;  
183.     }  
184.       
185.     public int clamp(int value) {  
186.         return value > 255 ? 255 :  
187.             (value < 0 ? 0 : value);  
188.     }  
189.       
190.     private int[] getPixel(double j, double k, int width, int height,  
191.             int[] inPixels) {  
192.         int row = (int) j;  
193.         int col = (int) k;  
194.         if (row >= height) {  
195.             row = height - 1;  
196.         }  
197.         if (row < 0) {  
198.             row = 0;  
199.         }  
200.         if (col < 0) {  
201.             col = 0;  
202.         }  
203.         if (col >= width) {  
204.             col = width - 1;  
205.         }  
206.         int index = row * width + col;  
207.         int[] rgb = new int[3];  
208.         rgb[0] = (inPixels[index] >> 16) & 0xff;  
209.         rgb[1] = (inPixels[index] >> 8) & 0xff;  
210.         rgb[2] = inPixels[index] & 0xff;  
211.         return rgb;  
212.     }  
213.     public BufferedImage createCompatibleDestImage(  
214.             BufferedImage src, ColorModel dstCM) {  
215.         if ( dstCM == null )  
216.             dstCM = src.getColorModel();  
217.         return new BufferedImage(dstCM,   
218.                 dstCM.createCompatibleWritableRaster(destW, destH),   
219.                 dstCM.isAlphaPremultiplied(), null);  
220.     }  
221. }  

运行效果：原图

双立方插值放大以后：

总结：

基于这里三种方法实现的双立方插值以后图片跟原图像相比，都有一定模糊

这里时候可以通过后续处理实现图像锐化与对比度提升即可得到Sharpen版本

当然也可以通过寻找更加合适的R(x)函数来实现双立方卷积插值过程时保留

图像边缘与对比度。

III.图像处理-OpenCV中resize函数五种插值算法的实现过程

（转载）

最新版OpenCV2.4.7中，cv::resize函数有五种插值算法：最近邻、双线性、双三次、基于像素区域关系、兰索斯插值。下面用for循环代替cv::resize函数来说明其详细的插值实现过程，其中部分代码摘自于cv::resize函数中的源代码。

每种插值算法的前部分代码是相同的，如下：

[cpp] view plaincopy

1. cv::Mat matSrc, matDst1, matDst2;  
2.   
3. matSrc = cv::imread("lena.jpg", 2 | 4);  
4. matDst1 = cv::Mat(cv::Size(800, 1000), matSrc.type(), cv::Scalar::all(0));  
5. matDst2 = cv::Mat(matDst1.size(), matSrc.type(), cv::Scalar::all(0));  
6.   
7. double scale_x = (double)matSrc.cols / matDst1.cols;  
8. double scale_y = (double)matSrc.rows / matDst1.rows;  

1、最近邻：公式，

                  

[cpp] view plaincopy

1. for (int i = 0; i < matDst1.cols; ++i)  
2. {  
3.     int sx = cvFloor(i * scale_x);  
4.     sx = std::min(sx, matSrc.cols - 1);  
5.     for (int j = 0; j < matDst1.rows; ++j)  
6.     {  
7.         int sy = cvFloor(j * scale_y);  
8.         sy = std::min(sy, matSrc.rows - 1);  
9.         matDst1.at<cv::Vec3b>(j, i) = matSrc.at<cv::Vec3b>(sy, sx);  
10.     }  
11. }  
12. cv::imwrite("nearest_1.jpg", matDst1);  
13.   
14. cv::resize(matSrc, matDst2, matDst1.size(), 0, 0, 0);  
15. cv::imwrite("nearest_2.jpg", matDst2);  

2、双线性：由相邻的四像素(2*2)计算得出，公式，

[cpp] view plaincopy

1. uchar* dataDst = matDst1.data;  
2. int stepDst = matDst1.step;  
3. uchar* dataSrc = matSrc.data;  
4. int stepSrc = matSrc.step;  
5. int iWidthSrc = matSrc.cols;  
6. int iHiehgtSrc = matSrc.rows;  
7.   
8. for (int j = 0; j < matDst1.rows; ++j)  
9. {  
10.     float fy = (float)((j + 0.5) * scale_y - 0.5);  
11.     int sy = cvFloor(fy);  
12.     fy -= sy;  
13.     sy = std::min(sy, iHiehgtSrc - 2);  
14.     sy = std::max(0, sy);  
15.   
16.     short cbufy[2];  
17.     cbufy[0] = cv::saturate_cast<short>((1.f - fy) * 2048);  
18.     cbufy[1] = 2048 - cbufy[0];  
19.   
20.     for (int i = 0; i < matDst1.cols; ++i)  
21.     {  
22.         float fx = (float)((i + 0.5) * scale_x - 0.5);  
23.         int sx = cvFloor(fx);  
24.         fx -= sx;  
25.   
26.         if (sx < 0) {  
27.             fx = 0, sx = 0;  
28.         }  
29.         if (sx >= iWidthSrc - 1) {  
30.             fx = 0, sx = iWidthSrc - 2;  
31.         }  
32.   
33.         short cbufx[2];  
34.         cbufx[0] = cv::saturate_cast<short>((1.f - fx) * 2048);  
35.         cbufx[1] = 2048 - cbufx[0];  
36.   
37.         for (int k = 0; k < matSrc.channels(); ++k)  
38.         {  
39.             *(dataDst+ j*stepDst + 3*i + k) = (*(dataSrc + sy*stepSrc + 3*sx + k) * cbufx[0] * cbufy[0] +   
40.                 *(dataSrc + (sy+1)*stepSrc + 3*sx + k) * cbufx[0] * cbufy[1] +   
41.                 *(dataSrc + sy*stepSrc + 3*(sx+1) + k) * cbufx[1] * cbufy[0] +   
42.                 *(dataSrc + (sy+1)*stepSrc + 3*(sx+1) + k) * cbufx[1] * cbufy[1]) >> 22;  
43.         }  
44.     }  
45. }  
46. cv::imwrite("linear_1.jpg", matDst1);  
47.   
48. cv::resize(matSrc, matDst2, matDst1.size(), 0, 0, 1);  
49. cv::imwrite("linear_2.jpg", matDst2);  

3、双三次：由相邻的4*4像素计算得出，公式类似于双线性

[cpp] view plaincopy

1. int iscale_x = cv::saturate_cast<int>(scale_x);  
2. int iscale_y = cv::saturate_cast<int>(scale_y);  
3.   
4. for (int j = 0; j < matDst1.rows; ++j)  
5. {  
6.     float fy = (float)((j + 0.5) * scale_y - 0.5);  
7.     int sy = cvFloor(fy);  
8.     fy -= sy;  
9.     sy = std::min(sy, matSrc.rows - 3);  
10.     sy = std::max(1, sy);  
11.   
12.     const float A = -0.75f;  
13.   
14.     float coeffsY[4];  
15.     coeffsY[0] = ((A*(fy + 1) - 5*A)*(fy + 1) + 8*A)*(fy + 1) - 4*A;  
16.     coeffsY[1] = ((A + 2)*fy - (A + 3))*fy*fy + 1;  
17.     coeffsY[2] = ((A + 2)*(1 - fy) - (A + 3))*(1 - fy)*(1 - fy) + 1;  
18.     coeffsY[3] = 1.f - coeffsY[0] - coeffsY[1] - coeffsY[2];  
19.   
20.     short cbufY[4];  
21.     cbufY[0] = cv::saturate_cast<short>(coeffsY[0] * 2048);  
22.     cbufY[1] = cv::saturate_cast<short>(coeffsY[1] * 2048);  
23.     cbufY[2] = cv::saturate_cast<short>(coeffsY[2] * 2048);  
24.     cbufY[3] = cv::saturate_cast<short>(coeffsY[3] * 2048);  
25.   
26.     for (int i = 0; i < matDst1.cols; ++i)  
27.     {  
28.         float fx = (float)((i + 0.5) * scale_x - 0.5);  
29.         int sx = cvFloor(fx);  
30.         fx -= sx;  
31.   
32.         if (sx < 1) {  
33.             fx = 0, sx = 1;  
34.         }  
35.         if (sx >= matSrc.cols - 3) {  
36.             fx = 0, sx = matSrc.cols - 3;  
37.         }  
38.   
39.         float coeffsX[4];  
40.         coeffsX[0] = ((A*(fx + 1) - 5*A)*(fx + 1) + 8*A)*(fx + 1) - 4*A;  
41.         coeffsX[1] = ((A + 2)*fx - (A + 3))*fx*fx + 1;  
42.         coeffsX[2] = ((A + 2)*(1 - fx) - (A + 3))*(1 - fx)*(1 - fx) + 1;  
43.         coeffsX[3] = 1.f - coeffsX[0] - coeffsX[1] - coeffsX[2];  
44.   
45.         short cbufX[4];  
46.         cbufX[0] = cv::saturate_cast<short>(coeffsX[0] * 2048);  
47.         cbufX[1] = cv::saturate_cast<short>(coeffsX[1] * 2048);  
48.         cbufX[2] = cv::saturate_cast<short>(coeffsX[2] * 2048);  
49.         cbufX[3] = cv::saturate_cast<short>(coeffsX[3] * 2048);  
50.   
51.         for (int k = 0; k < matSrc.channels(); ++k)  
52.         {  
53.             matDst1.at<cv::Vec3b>(j, i)[k] = abs((matSrc.at<cv::Vec3b>(sy-1, sx-1)[k] * cbufX[0] * cbufY[0] + matSrc.at<cv::Vec3b>(sy, sx-1)[k] * cbufX[0] * cbufY[1] +  
54.                 matSrc.at<cv::Vec3b>(sy+1, sx-1)[k] * cbufX[0] * cbufY[2] + matSrc.at<cv::Vec3b>(sy+2, sx-1)[k] * cbufX[0] * cbufY[3] +  
55.                 matSrc.at<cv::Vec3b>(sy-1, sx)[k] * cbufX[1] * cbufY[0] + matSrc.at<cv::Vec3b>(sy, sx)[k] * cbufX[1] * cbufY[1] +  
56.                 matSrc.at<cv::Vec3b>(sy+1, sx)[k] * cbufX[1] * cbufY[2] + matSrc.at<cv::Vec3b>(sy+2, sx)[k] * cbufX[1] * cbufY[3] +  
57.                 matSrc.at<cv::Vec3b>(sy-1, sx+1)[k] * cbufX[2] * cbufY[0] + matSrc.at<cv::Vec3b>(sy, sx+1)[k] * cbufX[2] * cbufY[1] +  
58.                 matSrc.at<cv::Vec3b>(sy+1, sx+1)[k] * cbufX[2] * cbufY[2] + matSrc.at<cv::Vec3b>(sy+2, sx+1)[k] * cbufX[2] * cbufY[3] +  
59.                 matSrc.at<cv::Vec3b>(sy-1, sx+2)[k] * cbufX[3] * cbufY[0] + matSrc.at<cv::Vec3b>(sy, sx+2)[k] * cbufX[3] * cbufY[1] +  
60.                 matSrc.at<cv::Vec3b>(sy+1, sx+2)[k] * cbufX[3] * cbufY[2] + matSrc.at<cv::Vec3b>(sy+2, sx+2)[k] * cbufX[3] * cbufY[3] ) >> 22);  
61.         }  
62.     }  
63. }  
64. cv::imwrite("cubic_1.jpg", matDst1);  
65.   
66. cv::resize(matSrc, matDst2, matDst1.size(), 0, 0, 2);  
67. cv::imwrite("cubic_2.jpg", matDst2);  

4、基于像素区域关系：共分三种情况，图像放大时类似于双线性插值，图像缩小(x轴、y轴同时缩小)又分两种情况，此情况下可以避免波纹出现。

[cpp] view plaincopy

1. cv::resize(matSrc, matDst2, matDst1.size(), 0, 0, 3);  
2. cv::imwrite("area_2.jpg", matDst2);  
3.   
4. double inv_scale_x = 1. / scale_x;  
5. double inv_scale_y = 1. / scale_y;  
6. int iscale_x = cv::saturate_cast<int>(scale_x);  
7. int iscale_y = cv::saturate_cast<int>(scale_y);  
8. bool is_area_fast = std::abs(scale_x - iscale_x) < DBL_EPSILON && std::abs(scale_y - iscale_y) < DBL_EPSILON;  
9.   
10. if (scale_x >= 1 && scale_y >= 1) //zoom out  
11. {  
12.     if (is_area_fast) //integer multiples  
13.     {  
14.         for (int j = 0; j < matDst1.rows; ++j)  
15.         {  
16.             int sy = j * scale_y;  
17.   
18.             for (int i = 0; i < matDst1.cols; ++i)  
19.             {  
20.                 int sx = i * scale_x;  
21.   
22.                 matDst1.at<cv::Vec3b>(j, i) = matSrc.at<cv::Vec3b>(sy, sx);  
23.             }  
24.         }  
25.         cv::imwrite("area_1.jpg", matDst1);  
26.         return 0;  
27.     }  
28.   
29.     for (int j = 0; j < matDst1.rows; ++j)  
30.     {  
31.         double fsy1 = j * scale_y;  
32.         double fsy2 = fsy1 + scale_y;  
33.         double cellHeight = cv::min(scale_y, matSrc.rows - fsy1);  
34.   
35.         int sy1 = cvCeil(fsy1), sy2 = cvFloor(fsy2);  
36.   
37.         sy2 = std::min(sy2, matSrc.rows - 1);  
38.         sy1 = std::min(sy1, sy2);  
39.   
40.         float cbufy[2];  
41.         cbufy[0] = (float)((sy1 - fsy1) / cellHeight);  
42.         cbufy[1] = (float)(std::min(std::min(fsy2 - sy2, 1.), cellHeight) / cellHeight);  
43.   
44.         for (int i = 0; i < matDst1.cols; ++i)  
45.         {  
46.             double fsx1 = i * scale_x;  
47.             double fsx2 = fsx1 + scale_x;  
48.             double cellWidth = std::min(scale_x, matSrc.cols - fsx1);  
49.   
50.             int sx1 = cvCeil(fsx1), sx2 = cvFloor(fsx2);  
51.   
52.             sx2 = std::min(sx2, matSrc.cols - 1);  
53.             sx1 = std::min(sx1, sx2);  
54.   
55.             float cbufx[2];  
56.             cbufx[0] = (float)((sx1 - fsx1) / cellWidth);  
57.             cbufx[1] = (float)(std::min(std::min(fsx2 - sx2, 1.), cellWidth) / cellWidth);  
58.   
59.             for (int k = 0; k < matSrc.channels(); ++k)  
60.             {  
61.                 matDst1.at<cv::Vec3b>(j, i)[k] = (uchar)(matSrc.at<cv::Vec3b>(sy1, sx1)[k] * cbufx[0] * cbufy[0] +   
62.                     matSrc.at<cv::Vec3b>(sy1 + 1, sx1)[k] * cbufx[0] * cbufy[1] +   
63.                     matSrc.at<cv::Vec3b>(sy1, sx1 + 1)[k] * cbufx[1] * cbufy[0] +   
64.                     matSrc.at<cv::Vec3b>(sy1 + 1, sx1 + 1)[k] * cbufx[1] * cbufy[1]);  
65.             }  
66.         }  
67.     }  
68.     cv::imwrite("area_1.jpg", matDst1);  
69.     return 0;  
70. }  
71.   
72. //zoom in,it is emulated using some variant of bilinear interpolation  
73. for (int j = 0; j < matDst1.rows; ++j)  
74. {  
75.     int  sy = cvFloor(j * scale_y);  
76.     float fy = (float)((j + 1) - (sy + 1) * inv_scale_y);  
77.     fy = fy <= 0 ? 0.f : fy - cvFloor(fy);  
78.   
79.     short cbufy[2];  
80.     cbufy[0] = cv::saturate_cast<short>((1.f - fy) * 2048);  
81.     cbufy[1] = 2048 - cbufy[0];  
82.   
83.     for (int i = 0; i < matDst1.cols; ++i)  
84.     {  
85.         int sx = cvFloor(i * scale_x);  
86.         float fx = (float)((i + 1) - (sx + 1) * inv_scale_x);  
87.         fx = fx < 0 ? 0.f : fx - cvFloor(fx);  
88.   
89.         if (sx < 0) {  
90.             fx = 0, sx = 0;  
91.         }  
92.   
93.         if (sx >= matSrc.cols - 1) {  
94.             fx = 0, sx = matSrc.cols - 2;  
95.         }  
96.   
97.         short cbufx[2];  
98.         cbufx[0] = cv::saturate_cast<short>((1.f - fx) * 2048);  
99.         cbufx[1] = 2048 - cbufx[0];  
100.   
101.         for (int k = 0; k < matSrc.channels(); ++k)  
102.         {  
103.             matDst1.at<cv::Vec3b>(j, i)[k] = (matSrc.at<cv::Vec3b>(sy, sx)[k] * cbufx[0] * cbufy[0] +   
104.                 matSrc.at<cv::Vec3b>(sy + 1, sx)[k] * cbufx[0] * cbufy[1] +   
105.                 matSrc.at<cv::Vec3b>(sy, sx + 1)[k] * cbufx[1] * cbufy[0] +   
106.                 matSrc.at<cv::Vec3b>(sy + 1, sx + 1)[k] * cbufx[1] * cbufy[1]) >> 22;  
107.         }  
108.     }  
109. }  
110. cv::imwrite("area_1.jpg", matDst1);  

5、兰索斯插值：由相邻的8*8像素计算得出，公式类似于双线性

[cpp] view plaincopy

1. int iscale_x = cv::saturate_cast<int>(scale_x);  
2. int iscale_y = cv::saturate_cast<int>(scale_y);  
3.   
4. for (int j = 0; j < matDst1.rows; ++j)  
5. {  
6.     float fy = (float)((j + 0.5) * scale_y - 0.5);  
7.     int sy = cvFloor(fy);  
8.     fy -= sy;  
9.     sy = std::min(sy, matSrc.rows - 5);  
10.     sy = std::max(3, sy);  
11.   
12.     const double s45 = 0.70710678118654752440084436210485;  
13.     const double cs[][2] = {{1, 0}, {-s45, -s45}, {0, 1}, {s45, -s45}, {-1, 0}, {s45, s45}, {0, -1}, {-s45, s45}};  
14.     float coeffsY[8];  
15.   
16.     if (fy < FLT_EPSILON) {  
17.         for (int t = 0; t < 8; t++)  
18.             coeffsY[t] = 0;  
19.         coeffsY[3] = 1;  
20.     } else {  
21.         float sum = 0;  
22.         double y0 = -(fy + 3) * CV_PI * 0.25, s0 = sin(y0), c0 = cos(y0);  
23.   
24.         for (int t = 0; t < 8; ++t)  
25.         {  
26.             double dy = -(fy + 3 -t) * CV_PI * 0.25;  
27.             coeffsY[t] = (float)((cs[t][0] * s0 + cs[t][1] * c0) / (dy * dy));  
28.             sum += coeffsY[t];  
29.         }  
30.   
31.         sum = 1.f / sum;  
32.         for (int t = 0; t < 8; ++t)  
33.             coeffsY[t] *= sum;  
34.     }  
35.   
36.     short cbufY[8];  
37.     cbufY[0] = cv::saturate_cast<short>(coeffsY[0] * 2048);  
38.     cbufY[1] = cv::saturate_cast<short>(coeffsY[1] * 2048);  
39.     cbufY[2] = cv::saturate_cast<short>(coeffsY[2] * 2048);  
40.     cbufY[3] = cv::saturate_cast<short>(coeffsY[3] * 2048);  
41.     cbufY[4] = cv::saturate_cast<short>(coeffsY[4] * 2048);  
42.     cbufY[5] = cv::saturate_cast<short>(coeffsY[5] * 2048);  
43.     cbufY[6] = cv::saturate_cast<short>(coeffsY[6] * 2048);  
44.     cbufY[7] = cv::saturate_cast<short>(coeffsY[7] * 2048);  
45.   
46.     for (int i = 0; i < matDst1.cols; ++i)  
47.     {  
48.         float fx = (float)((i + 0.5) * scale_x - 0.5);  
49.         int sx = cvFloor(fx);  
50.         fx -= sx;  
51.   
52.         if (sx < 3) {  
53.             fx = 0, sx = 3;  
54.         }  
55.         if (sx >= matSrc.cols - 5) {  
56.             fx = 0, sx = matSrc.cols - 5;  
57.         }  
58.   
59.         float coeffsX[8];  
60.   
61.         if (fx < FLT_EPSILON) {  
62.             for ( int t = 0; t < 8; t++ )  
63.                 coeffsX[t] = 0;  
64.             coeffsX[3] = 1;  
65.         } else {  
66.             float sum = 0;  
67.             double x0 = -(fx + 3) * CV_PI * 0.25, s0 = sin(x0), c0 = cos(x0);  
68.   
69.             for (int t = 0; t < 8; ++t)  
70.             {  
71.                 double dx = -(fx + 3 -t) * CV_PI * 0.25;  
72.                 coeffsX[t] = (float)((cs[t][0] * s0 + cs[t][1] * c0) / (dx * dx));  
73.                 sum += coeffsX[t];  
74.             }  
75.   
76.             sum = 1.f / sum;  
77.             for (int t = 0; t < 8; ++t)  
78.                 coeffsX[t] *= sum;  
79.         }  
80.   
81.         short cbufX[8];  
82.         cbufX[0] = cv::saturate_cast<short>(coeffsX[0] * 2048);  
83.         cbufX[1] = cv::saturate_cast<short>(coeffsX[1] * 2048);  
84.         cbufX[2] = cv::saturate_cast<short>(coeffsX[2] * 2048);  
85.         cbufX[3] = cv::saturate_cast<short>(coeffsX[3] * 2048);  
86.         cbufX[4] = cv::saturate_cast<short>(coeffsX[4] * 2048);  
87.         cbufX[5] = cv::saturate_cast<short>(coeffsX[5] * 2048);  
88.         cbufX[6] = cv::saturate_cast<short>(coeffsX[6] * 2048);  
89.         cbufX[7] = cv::saturate_cast<short>(coeffsX[7] * 2048);  
90.   
91.         for (int k = 0; k < matSrc.channels(); ++k)  
92.         {  
93.             matDst1.at<cv::Vec3b>(j, i)[k] = abs((matSrc.at<cv::Vec3b>(sy-3, sx-3)[k] * cbufX[0] * cbufY[0] + matSrc.at<cv::Vec3b>(sy-2, sx-3)[k] * cbufX[0] * cbufY[1] +  
94.                 matSrc.at<cv::Vec3b>(sy-1, sx-3)[k] * cbufX[0] * cbufY[2] + matSrc.at<cv::Vec3b>(sy, sx-3)[k] * cbufX[0] * cbufY[3] +  
95.                 matSrc.at<cv::Vec3b>(sy+1, sx-3)[k] * cbufX[0] * cbufY[4] + matSrc.at<cv::Vec3b>(sy+2, sx-3)[k] * cbufX[0] * cbufY[5] +  
96.                 matSrc.at<cv::Vec3b>(sy+3, sx-3)[k] * cbufX[0] * cbufY[6] + matSrc.at<cv::Vec3b>(sy+4, sx-3)[k] * cbufX[0] * cbufY[7] +  
97.   
98.                 matSrc.at<cv::Vec3b>(sy-3, sx-2)[k] * cbufX[1] * cbufY[0] + matSrc.at<cv::Vec3b>(sy-2, sx-2)[k] * cbufX[1] * cbufY[1] +  
99.                 matSrc.at<cv::Vec3b>(sy-1, sx-2)[k] * cbufX[1] * cbufY[2] + matSrc.at<cv::Vec3b>(sy, sx-2)[k] * cbufX[1] * cbufY[3] +  
100.                 matSrc.at<cv::Vec3b>(sy+1, sx-2)[k] * cbufX[1] * cbufY[4] + matSrc.at<cv::Vec3b>(sy+2, sx-2)[k] * cbufX[1] * cbufY[5] +  
101.                 matSrc.at<cv::Vec3b>(sy+3, sx-2)[k] * cbufX[1] * cbufY[6] + matSrc.at<cv::Vec3b>(sy+4, sx-2)[k] * cbufX[1] * cbufY[7] +  
102.   
103.                 matSrc.at<cv::Vec3b>(sy-3, sx-1)[k] * cbufX[2] * cbufY[0] + matSrc.at<cv::Vec3b>(sy-2, sx-1)[k] * cbufX[2] * cbufY[1] +  
104.                 matSrc.at<cv::Vec3b>(sy-1, sx-1)[k] * cbufX[2] * cbufY[2] + matSrc.at<cv::Vec3b>(sy, sx-1)[k] * cbufX[2] * cbufY[3] +  
105.                 matSrc.at<cv::Vec3b>(sy+1, sx-1)[k] * cbufX[2] * cbufY[4] + matSrc.at<cv::Vec3b>(sy+2, sx-1)[k] * cbufX[2] * cbufY[5] +  
106.                 matSrc.at<cv::Vec3b>(sy+3, sx-1)[k] * cbufX[2] * cbufY[6] + matSrc.at<cv::Vec3b>(sy+4, sx-1)[k] * cbufX[2] * cbufY[7] +  
107.   
108.                 matSrc.at<cv::Vec3b>(sy-3, sx)[k] * cbufX[3] * cbufY[0] + matSrc.at<cv::Vec3b>(sy-2, sx)[k] * cbufX[3] * cbufY[1] +  
109.                 matSrc.at<cv::Vec3b>(sy-1, sx)[k] * cbufX[3] * cbufY[2] + matSrc.at<cv::Vec3b>(sy, sx)[k] * cbufX[3] * cbufY[3] +  
110.                 matSrc.at<cv::Vec3b>(sy+1, sx)[k] * cbufX[3] * cbufY[4] + matSrc.at<cv::Vec3b>(sy+2, sx)[k] * cbufX[3] * cbufY[5] +  
111.                 matSrc.at<cv::Vec3b>(sy+3, sx)[k] * cbufX[3] * cbufY[6] + matSrc.at<cv::Vec3b>(sy+4, sx)[k] * cbufX[3] * cbufY[7] +  
112.   
113.                 matSrc.at<cv::Vec3b>(sy-3, sx+1)[k] * cbufX[4] * cbufY[0] + matSrc.at<cv::Vec3b>(sy-2, sx+1)[k] * cbufX[4] * cbufY[1] +  
114.                 matSrc.at<cv::Vec3b>(sy-1, sx+1)[k] * cbufX[4] * cbufY[2] + matSrc.at<cv::Vec3b>(sy, sx+1)[k] * cbufX[4] * cbufY[3] +  
115.                 matSrc.at<cv::Vec3b>(sy+1, sx+1)[k] * cbufX[4] * cbufY[4] + matSrc.at<cv::Vec3b>(sy+2, sx+1)[k] * cbufX[4] * cbufY[5] +  
116.                 matSrc.at<cv::Vec3b>(sy+3, sx+1)[k] * cbufX[4] * cbufY[6] + matSrc.at<cv::Vec3b>(sy+4, sx+1)[k] * cbufX[4] * cbufY[7] +  
117.   
118.                 matSrc.at<cv::Vec3b>(sy-3, sx+2)[k] * cbufX[5] * cbufY[0] + matSrc.at<cv::Vec3b>(sy-2, sx+2)[k] * cbufX[5] * cbufY[1] +  
119.                 matSrc.at<cv::Vec3b>(sy-1, sx+2)[k] * cbufX[5] * cbufY[2] + matSrc.at<cv::Vec3b>(sy, sx+2)[k] * cbufX[5] * cbufY[3] +  
120.                 matSrc.at<cv::Vec3b>(sy+1, sx+2)[k] * cbufX[5] * cbufY[4] + matSrc.at<cv::Vec3b>(sy+2, sx+2)[k] * cbufX[5] * cbufY[5] +  
121.                 matSrc.at<cv::Vec3b>(sy+3, sx+2)[k] * cbufX[5] * cbufY[6] + matSrc.at<cv::Vec3b>(sy+4, sx+2)[k] * cbufX[5] * cbufY[7] +  
122.   
123.                 matSrc.at<cv::Vec3b>(sy-3, sx+3)[k] * cbufX[6] * cbufY[0] + matSrc.at<cv::Vec3b>(sy-2, sx+3)[k] * cbufX[6] * cbufY[1] +  
124.                 matSrc.at<cv::Vec3b>(sy-1, sx+3)[k] * cbufX[6] * cbufY[2] + matSrc.at<cv::Vec3b>(sy, sx+3)[k] * cbufX[6] * cbufY[3] +  
125.                 matSrc.at<cv::Vec3b>(sy+1, sx+3)[k] * cbufX[6] * cbufY[4] + matSrc.at<cv::Vec3b>(sy+2, sx+3)[k] * cbufX[6] * cbufY[5] +  
126.                 matSrc.at<cv::Vec3b>(sy+3, sx+3)[k] * cbufX[6] * cbufY[6] + matSrc.at<cv::Vec3b>(sy+4, sx+3)[k] * cbufX[6] * cbufY[7] +  
127.   
128.                 matSrc.at<cv::Vec3b>(sy-3, sx+4)[k] * cbufX[7] * cbufY[0] + matSrc.at<cv::Vec3b>(sy-2, sx+4)[k] * cbufX[7] * cbufY[1] +  
129.                 matSrc.at<cv::Vec3b>(sy-1, sx+4)[k] * cbufX[7] * cbufY[2] + matSrc.at<cv::Vec3b>(sy, sx+4)[k] * cbufX[7] * cbufY[3] +  
130.                 matSrc.at<cv::Vec3b>(sy+1, sx+4)[k] * cbufX[7] * cbufY[4] + matSrc.at<cv::Vec3b>(sy+2, sx+4)[k] * cbufX[7] * cbufY[5] +  
131.                 matSrc.at<cv::Vec3b>(sy+3, sx+4)[k] * cbufX[7] * cbufY[6] + matSrc.at<cv::Vec3b>(sy+4, sx+4)[k] * cbufX[7] * cbufY[7] ) >> 22);// 4194304  
132.         }  
133.     }  
134. }  
135. cv::imwrite("Lanczos_1.jpg", matDst1);  
136.   
137. cv::resize(matSrc, matDst2, matDst1.size(), 0, 0, 4);  
138. cv::imwrite("Lanczos_2.jpg", matDst2);  

    以上代码的实现结果与cv::resize函数相同，但是执行效率非常低，只是为了详细说明插值过程。OpenCV中默认采用C++ Concurrency进行优化加速，你也可以采用TBB、OpenMP等进行优化加速。