怎么将一个矩形图像变形为一个梯形( 或者说是任意四边形)？有源代码的话最好

如何将一个矩形图像变形为一个梯形( 或者说是任意四边形)？有源代码的话最好。
如何将一个矩形图像变形为一个梯形( 或者说是任意四边形)？有源代码的话最好。




------解决方案--------------------
// 白像素
#define HIGH 255

// 黑像素
#define LOW 0

// π值
#define PI 3.14159265f

/*------matrix--------矩阵计算------------------------------
l: 输入矩阵1
m: 输入矩阵2
n: 输出矩阵
-----------------------*/
inline void matrix(double l[4][4], double m[4][4], double n[4][4])
{
int i, j, k;
double p;

for(i = 0; i < 4; i++)
{
for(j = 0; j < 4; j++)
{
p = 0;
for(k = 0; k < 4; k++)
p += l[i][k] * m[k][j];
n[i][j] = p;
}
}
}

/*------param_pers-------计算透视变换的参数-----------------
xsize: 图像宽度
ysize: 图像高度
k: 变换参数
a: 放大率（x方向）
b: 放大率（y方向）
x0: 移动量（x方向）
y0: 移动量（y方向）
z0: 移动量（z方向）
x: 回转角（x方向）
y: 回转角（y方向）
z: 回转角（z方向）
t: 视点位置（z）
s: 屏幕位置（z）
-----------------------*/
inline void param_pers(int xsize, int ysize, float k[9], float a, float b, float x0, float y0, float z0, float x, float y, float z, float t, float s)
{
double l[4][4], m[4][4], n[4][4];
float k1, k2, k3, k4, k5, k6, k7, k8, k9;
double u, v, w;
int xs = xsize / 2;
int ys = ysize / 2;

u = x * PI / 180.0; v = y * PI / 180.0; w = z * PI / 180.0;
l[0][0] = 1.0 / xs; l[0][1] = 0 ; l[0][2] = 0 ; l[0][3] = 0 ;
l[1][0] = 0 ; l[1][1] = -1.0 / xs; l[1][2] = 0 ; l[1][3] = 0 ;
l[2][0] = 0 ; l[2][1] = 0 ; l[2][2] = 1 ; l[2][3] = 0 ;
l[3][0] = 0 ; l[3][1] = 0 ; l[3][2] = 0 ; l[3][3] = 1 ;
m[0][0] = a ; m[0][1] = 0 ; m[0][2] = 0 ; m[0][3] = 0 ;
m[1][0] = 0 ; m[1][1] = b ; m[1][2] = 0 ; m[1][3] = 0 ;
m[2][0] = 0 ; m[2][1] = 0 ; m[2][2] = 1 ; m[2][3] = 0 ;
m[3][0] = 0 ; m[3][1] = 0 ; m[3][2] = 0 ; m[3][3] = 1 ;
matrix(l, m, n);

l[0][0] = 1 ; l[0][1] = 0 ; l[0][2] = 0 ; l[0][3] = 0 ;
l[1][0] = 0 ; l[1][1] = 1 ; l[1][2] = 0 ; l[1][3] = 0 ;
l[2][0] = 0 ; l[2][1] = 0 ; l[2][2] = 1 ; l[2][3] = 0 ;
l[3][0] = x0 ; l[3][1] = y0 ; l[3][2] = z0 ; l[3][3] = 1 ;
matrix(n, l, m);

n[0][0] = cos(w) ; n[0][1] = sin(w) ; n[0][2] = 0 ; n[0][3] = 0 ;
n[1][0] = -sin(w) ; n[1][1] = cos(w) ; n[1][2] = 0 ; n[1][3] = 0 ;
n[2][0] = 0 ; n[2][1] = 0 ; n[2][2] = 1 ; n[2][3] = 0 ;
n[3][0] = 0 ; n[3][1] = 0 ; n[3][2] = 0 ; n[3][3] = 1 ;
matrix(m, n, l);

m[0][0] = 1 ; m[0][1] = 0 ; m[0][2] = 0 ; m[0][3] = 0 ;
m[1][0] = 0 ; m[1][1] = cos(u) ; m[1][2] = sin(u) ; m[1][3] = 0 ;
m[2][0] = 0 ; m[2][1] = -sin(u) ; m[2][2] = cos(u) ; m[2][3] = 0 ;
m[3][0] = 0 ; m[3][1] = 0 ; m[3][2] = 0 ; m[3][3] = 1 ;
matrix(l, m, n);

l[0][0] = cos(v) ; l[0][1] = 0 ; l[0][2] = sin(v) ; l[0][3] = 0 ;
l[1][0] = 0 ; l[1][1] = 1 ; l[1][2] = 0 ; l[1][3] = 0 ;
l[2][0] = -sin(v) ; l[2][1] = 0 ; l[2][2] = cos(v) ; l[2][3] = 0 ;
l[3][0] = 0 ; l[3][1] = 0 ; l[3][2] = 0 ; l[3][3] = 1 ;
matrix(n, l, m);

n[0][0] = 1 ; n[0][1] = 0 ; n[0][2] = 0 ; n[0][3] = 0 ;
n[1][0] = 0 ; n[1][1] = 1 ; n[1][2] = 0 ; n[1][3] = 0 ;
n[2][0] = 0 ; n[2][1] = 0 ; n[2][2] = -1 ; n[2][3] = 0 ;