双线插值算法原理
(原理部分转自:http://blog.163.com/guohuanhuan_cool@126/blog/static/167614238201161525538402/)
图 像的缩放很好理解,就是图像的放大和缩小。传统的绘画工具中,有一种叫做“放大尺”的绘画工具,画家常用它来放大图画。当然,在计算机上,我们不再需要用 放大尺去放大或缩小图像了,把这个工作交给程序来完成就可以了。下面就来讲讲计算机怎么来放大缩小图象;在本文中,我们所说的图像都是指点阵图,也就是用 一个像素矩阵来描述图像的方法,对于另一种图像:用函数来描述图像的矢量图,不在本文讨论之列。
越是简单的模型越适合用来举例子,我们就举个简单 的图像: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也体现出了这一特点;
双线插值算法实现(C/C++)
typedef long fixedpoint_t;/* fixed point type,低8位定义为小数部分 */
#define itofx(x) ((x) << 8) /* Integer to fixed point */
#define ftofx(x) (long)((x) * 256) /* Float to fixed point */
#define dtofx(x) (long)((x) * 256) /* Double to fixed point */
#define fxtoi(x) ((x) >> 8) /* Fixed point to integer */
#define fxtof(x) ((float)(x) / 256) /* Fixed point to float */
#define fxtod(x) ((double)(x) / 256) /* Fixed point to double */
#define mulfx(x,y) (((x) * (y)) >> 8) /* Multiply a fixed by a fixed */
#define divfx(x,y) (((x) << 8) / (y)) /* Divide a fixed by a fixed */
#define MAX(a,b) ((a)<(b)?(b):(a))
#ifndef NULL
#define NULL 0
#endif
#ifdef _WIN32
#define inline __inline
#endif
static inline unsigned short getPixelOverflow(const unsigned short * tpixels,const int pitch,const int sx,const int sy,const int x,const int y)
{
if(y<0)
return tpixels[x];
if(y>(sy-1))
return tpixels[(sy-1)*pitch+x];
if(x<0)
return tpixels[y*pitch];
if(x>(sx-1))
return tpixels[y*pitch+sx-1];
return tpixels[y*pitch+x];
}
/* only scale rgb555 to rgb555 */
unsigned short* imageScaleBilinear(const unsigned short*tpixels, const unsigned int sx,const unsigned int sy,const unsigned int new_width, const unsigned int new_height)
{
long dst_w = MAX(1, new_width);
long dst_h = MAX(1, new_height);
float dx = (float)sx / (float)dst_w;
float dy = (float)sy / (float)dst_h;
fixedpoint_t f_dx = ftofx(dx);
fixedpoint_t f_dy = ftofx(dy);
fixedpoint_t f_1 = itofx(1);
unsigned int pitch = ((sx*2+3)&0xfffffffc)/2;//tpixels type is unsigned short
/*
if source image data format is not rgb555,need recalculate "pitch" value.
eg:
rgb8888:pitch = (w * 4 + 3) & 0xfffffffc;if tpixels type is uint32_t,pitch/=4;
rgb888: pitch = (w * 3 + 3) & 0xfffffffc;
rgb8: pitch = (w * 1 + 3) & 0xfffffffc;
*/
int dst_offset_h;
int dst_offset_v = 0;
long i;
unsigned int new_pitch = (new_width*2+3)&0xfffffffc;
unsigned short* new_img;
new_img = new unsigned short[ new_pitch * new_height ];
if (!new_img){
return NULL;
}
new_pitch /= 2;
for (i=0; i < dst_h; i++)
{
long j;
dst_offset_h = 0;
for (j=0; j < dst_w; j++) {
/* Update bitmap */
fixedpoint_t f_i = itofx(i);
fixedpoint_t f_j = itofx(j);
fixedpoint_t f_a = mulfx(f_i, f_dy);
fixedpoint_t f_b = mulfx(f_j, f_dx);
const long m = fxtoi(f_a);
const long n = fxtoi(f_b);
fixedpoint_t f_f = f_a - itofx(m);
fixedpoint_t f_g = f_b - itofx(n);
const fixedpoint_t f_w1 = mulfx(f_1-f_f, f_1-f_g);
const fixedpoint_t f_w2 = mulfx(f_1-f_f, f_g);
const fixedpoint_t f_w3 = mulfx(f_f, f_1-f_g);
const fixedpoint_t f_w4 = mulfx(f_f, f_g);
unsigned short pixel1;
unsigned short pixel2;
unsigned short pixel3;
unsigned short pixel4;
register fixedpoint_t f_r1, f_r2, f_r3, f_r4,
f_g1, f_g2, f_g3, f_g4,
f_b1, f_b2, f_b3, f_b4;
/* nothing gets outside anyway */
pixel1 = getPixelOverflow(tpixels, pitch, sx, sy, n, m);
pixel2 = getPixelOverflow(tpixels, pitch, sx, sy, n + 1, m);
pixel3 = getPixelOverflow(tpixels, pitch, sx, sy, n, m + 1);
pixel4 = getPixelOverflow(tpixels, pitch, sx, sy, n + 1, m + 1);
#define RGB555GetRed(c) (((c) & 0x7C00) >> 10)
#define RGB555GetGreen(c) (((c) & 0x03E0) >> 5)
#define RGB555GetBlue(c) (((c) & 0x001F) >> 0)
f_r1 = itofx(RGB555GetRed(pixel1));
f_r2 = itofx(RGB555GetRed(pixel2));
f_r3 = itofx(RGB555GetRed(pixel3));
f_r4 = itofx(RGB555GetRed(pixel4));
f_g1 = itofx(RGB555GetGreen(pixel1));
f_g2 = itofx(RGB555GetGreen(pixel2));
f_g3 = itofx(RGB555GetGreen(pixel3));
f_g4 = itofx(RGB555GetGreen(pixel4));
f_b1 = itofx(RGB555GetBlue(pixel1));
f_b2 = itofx(RGB555GetBlue(pixel2));
f_b3 = itofx(RGB555GetBlue(pixel3));
f_b4 = itofx(RGB555GetBlue(pixel4));
#undef RGB555GetRed
#undef RGB555GetGreen
#undef RGB555GetBlue
{
const char red = (char) fxtoi(mulfx(f_w1, f_r1) + mulfx(f_w2, f_r2) + mulfx(f_w3, f_r3) + mulfx(f_w4, f_r4));
const char green = (char) fxtoi(mulfx(f_w1, f_g1) + mulfx(f_w2, f_g2) + mulfx(f_w3, f_g3) + mulfx(f_w4, f_g4));
const char blue = (char) fxtoi(mulfx(f_w1, f_b1) + mulfx(f_w2, f_b2) + mulfx(f_w3, f_b3) + mulfx(f_w4, f_b4));
#define RGB555(r,g,b) (((r) << 10) | ((g) << 5) | (b))
new_img[new_pitch*dst_offset_v+dst_offset_h] = RGB555(red,green,blue);
#undef RGB555
}
dst_offset_h++;
}
dst_offset_v++;
}
return new_img;
}
/**********************for test***************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#ifndef _WIN32
typedef long BOOL;
typedef long LONG;
typedef unsigned char BYTE;
typedef unsigned long DWORD;
typedef unsigned short WORD;
typedef struct {
WORD bfType;//2
DWORD bfSize;//4
WORD bfReserved1;//2
WORD bfReserved2;//2
DWORD bfOffBits;//4
}__attribute__((packed))BmpFileHead;
typedef struct{
DWORD biSize;//4
LONG biWidth;//4
LONG biHeight;//4
WORD biPlanes;//2
WORD biBitCount;//2
DWORD biCompression;//4
DWORD biSizeImage;//4
LONG biXPelsPerMeter;//4
LONG biYPelsPerMeter;//4
DWORD biClrUsed;//4
DWORD biClrImportant;//4
}__attribute__((packed))BmpInfohead;
#else
#include <windows.h>
typedef BITMAPFILEHEADER BmpFileHead;
typedef BITMAPINFOHEADER BmpInfohead;
#endif
/* only get 24bit bmp data,no bmp header*/
unsigned char * get_bmp24_data(const char * src_bmp24_file,int * sx,int * sy)
{
FILE *fp = fopen(src_bmp24_file, "rb");//
unsigned char * bmp_data = NULL;
if(fp)
{
BmpFileHead bmp_head;
BmpInfohead bmp_info;
unsigned int w,h,data_size;
fread(&bmp_head,sizeof(BmpFileHead),1,fp);
fread(&bmp_info,sizeof(BmpInfohead),1,fp);
w = bmp_info.biWidth;
h = bmp_info.biHeight;
data_size = bmp_head.bfSize - bmp_head.bfOffBits;
if(h<0)h=-h;
if(data_size != ((w*3+3)&0xfffffffc) * h )
{
printf("not 24 bit bmp file,file size = %u,w=%u,h=%u\n",data_size,w,h);
fclose(fp);
return NULL;
}
bmp_data = new unsigned char[data_size];
memset(bmp_data,0,data_size);
if(bmp_data)
{
fread(bmp_data,1,data_size,fp);
}
fclose(fp);
fp = NULL;
*sx = w;
*sy = h;
}
return bmp_data;
}
/* only save 24bit bmp file */
int save_bmp24(const char * filename,int width,int height,unsigned char *data)
{
BmpFileHead bmp_head;
BmpInfohead bmp_info;
int size = ((width*3+3)&0xfffffffc)*height;
FILE *fp = fopen(filename,"wb");
if(!fp)
{
perror("open file error");
return -1;
}
bmp_head.bfType=0x4d42;
bmp_head.bfSize=size+sizeof(BmpFileHead)+sizeof(BmpInfohead);//24+head+info no quad
bmp_head.bfReserved1=bmp_head.bfReserved2=0;
bmp_head.bfOffBits=bmp_head.bfSize-size;
bmp_info.biSize=40;
bmp_info.biWidth=width;
bmp_info.biHeight=height;
bmp_info.biPlanes=1;
bmp_info.biBitCount = 24;
bmp_info.biCompression=0;
bmp_info.biSizeImage=size;
bmp_info.biXPelsPerMeter=0;
bmp_info.biYPelsPerMeter=0;
bmp_info.biClrUsed=0;
bmp_info.biClrImportant=0;
fwrite(&bmp_head,1,sizeof(BmpFileHead),fp);
fwrite(&bmp_info,1,sizeof(BmpInfohead),fp);
fwrite(data,1,size,fp);
fclose(fp);
return 0;
}
/* convert rgb888 to rgb555 (RGB order:R G B) */
unsigned short * rgb888_rgb555(const unsigned char * rgb888,unsigned int width,unsigned int height)
{
unsigned int pitch = (width*2+3)&0xfffffffc;
unsigned int rgb888pitch = (width*3+3)&0xfffffffc;
unsigned short * rgb555=new unsigned short[pitch * height ];
if(rgb555)
{
unsigned int w,h;
pitch/=2;
for(h=0;h<height;h++)
{
unsigned char *p=(unsigned char *)&rgb888[rgb888pitch*h];
for(w=0;w<width;w++)
{
#define RGB555(r,g,b) ( (((r)>>3)<<10)|(((g)>>3)<<5)|(((b)>>3)) )
rgb555[h*pitch+w]=RGB555(*p,*(p+1),*(p+2));
#undef RGB555
p+=3;
}
}
}
return rgb555;
}
/* convert rgb555 to rgb888 (RGB order:R G B) */
unsigned char * rgb555_rgb888(const unsigned short * rgb555,unsigned int width,unsigned int height)
{
unsigned int pitch = (width*3+3)&0xfffffffc;
unsigned int rgb555pitch=(width*2+3)&0xfffffffc;
unsigned char * rgb888=new unsigned char[pitch * height];
if(rgb888)
{
#define RGB555GetRed(c) (((c) & 0x7C00) >> 10)
#define RGB555GetGreen(c) (((c) & 0x03E0) >> 5)
#define RGB555GetBlue(c) (((c) & 0x001F) )
unsigned int w,h;
rgb555pitch/=2;//rgb555为unsigned short
for(h=0;h<height;h++)
{
unsigned char * p = &rgb888[pitch*h];
for(w=0;w<width;w++)
{
unsigned short c = rgb555[h*rgb555pitch+w];
*p++=RGB555GetRed(c)<<3;
*p++=RGB555GetGreen(c)<<3;
*p++=RGB555GetBlue(c)<<3;
}
}
#undef RGB555GetRed
#undef RGB555GetGreen
#undef RGB555GetBlue
}
return rgb888;
}
int main(int argc,char * argv[])
{
/*
load rgb888 from file
---> convert rgb888 to rgb555
---> scale rgb555 to rgb555
---> conver rgb555 to rgb888
---> save to file
*/
int src_sx=0,src_sy=0;
unsigned char * src_bmp;
if(argc != 5)
{
USAGE:
printf("USAGE:%s srcimg dstimg dstimg_width dstimg_height\n",argv[0]);
return -1;
}
if(atoi(argv[3])<=0 || atoi(argv[4])<= 0)
{
goto USAGE;
}
src_bmp = get_bmp24_data(argv[1],&src_sx,&src_sy);
if(src_bmp)
{
unsigned short * src_rgb555=NULL;
printf("Load source file success[%dx%d]...\n",src_sx,src_sy);
src_rgb555 = rgb888_rgb555(src_bmp,src_sx,src_sy);
if(src_rgb555)
{
unsigned short * dst_rgb555;
int dst_w=atoi(argv[3]),dst_h=atoi(argv[4]);
struct timespec tpstart,tpend;
double timedif;
printf("Conver to RGB555 success...\n");
clock_gettime(CLOCK_MONOTONIC, &tpstart);
dst_rgb555 = imageScaleBilinear(src_rgb555,src_sx,src_sy,dst_w,dst_h);
clock_gettime(CLOCK_MONOTONIC, &tpend);
timedif = 1000.0L*(double)(tpend.tv_sec-tpstart.tv_sec)+(double)(tpend.tv_nsec-tpstart.tv_nsec)/1000000.0L;
fprintf(stderr,"Exec time:%.3lfms,\n",timedif);
if(dst_rgb555)
{
unsigned char * dst_rgb888;
printf("Scale to %dx%d success...\n",dst_w,dst_h);
dst_rgb888 = rgb555_rgb888(dst_rgb555,dst_w,dst_h);
if(dst_rgb888)
{
printf("Conver to RGB888 success...\n");
if(0==save_bmp24(argv[2],dst_w,dst_h,dst_rgb888))
{
printf("save to %s success...\n",argv[2]);
}
delete []dst_rgb888;
}
delete []dst_rgb555;
}
delete []src_rgb555;
}
delete []src_bmp;
}
return 0;
}
/***
*test result(ARM CPU,400MHz):
*src image size:1280x800
*dst image size:128x128,exec time:11 ~ 13 ms
*dst image size:80x80,exec time:5 ~ 6 ms
*dst image size:64x64,exec time:3 ~ 4 ms
*dst image size:48x48,exec time:2 ~ 3 ms
*dst image size:40x40,exec time:1 ~ 2 ms
*dst image size:32x32,exec time:about 1 ms
*/
/**********************test end***************************/