1透视投影变换概述
透视投影变换(Perspective Transformation)是对图像做视平面(Viewing Plane)的变换,通俗来讲可以用下面一组图来表示
图1 图中画了一只小鸟v
在图中有一幅画,画中是一只小鸟,小蓝人在画的左边,离画较远;小红人站在画的右边,离画较近。他们眼中的图分别是
图2 小蓝人眼中的画(左) 小红人眼中的画(右)
如果小蓝人想知道小红人看到的画是什么样的,那么他有两种方法,一种是他走到小红人站的位置,也就是说在小红人所处的视觉坐标下观察;另个一种方法就是将小蓝人看到的图像所处的视平面变换到小红人所处的视平面下。对于墙壁上的画来说,这两种方法的结果是一样的。如果换做是雕像那种三维的物体,结果就不一样了,但在精度要求不高的情况下,也是可以近似的v。
2透视投影变换原理
具体详细的原理的介绍,相关书籍中写的很详细,这里简单介绍下:
透视变换将(x,y)点映射为(u,v)点,如下所示:
u=(ax+by+c)/(dx+ey+f)
v=(gx+hy+1)/(dx+ey+f)
在物理上我们将2D图像嵌入到3D空间中,在包含2D图像的平面上应用透视变换,然后应用透视图回顾一下仿射变换的形式:
u=ax+by+c
v=dx+ey+f
u的等式中的参数a,b,c不出现在v的等式中,但是对于透视变换,参数d,e,f出现在两个方程中,u,v的分母相同,dx + ey + f是将3D透视变换应用于包含2D图像的平面之后的点的深度。在物理上我们增加了深度的维度。
假设我们有四个控制点p1到p4分别映射到矩形的四个顶点q1到q4(也可以不是矩形,任意四边形都是可以的)。p1到p4为原图中期望矫正成为矩形的四边形的顶点。这可以以矩阵形式写成
Ax=B
其中x是未知数的向量
x=(a,b,c,d,e,f,g,h)^T
且A是由下式给出 (文章使用word拟的草稿,公式不能粘贴,只能截图了)
B为
B=(q1.x q1.y q2.x q2.y q3.x q3.y q4.x q4.y)^T
解以上方程可以得到参数a,b,c,d,e,f,g,h。
利用变换阵
对图像进行变换即可。
3透视投影变换的实现
直接上代码
先是MATLAB版的,用于验证算法。后面还有C语言版的,比较实用。
MATLAB版:
img= imread('20160303.bmp');%读取图像,灰度
%img= rgb2gray(img);%如果是彩色的,要转成灰度!!!
[row col] = size(img);%输入图像的尺寸
%按照[行,列],而不是[x,y]
% 0----------------------->col
% |
% |
% |
% |
% |
% V
% row
p1=[21 57];%这四个点是根据本程序中使用的图像事先取好的,可以是图像中的任意的部分。顺序为左上,右上,左下,右下
p2=[21 108];
p3=[53 39];
p4=[53 120];
%得到p1,p2,p3,p4这个小矩形所对应的变换后的矩形的宽
w=round(sqrt((p1(1)-p2(1))^2+(p1(2)-p2(2))^2));
%矩形的高
h=round(sqrt((p1(1)-p3(1))^2+(p1(2)-p3(2))^2));
%这里是新的顶点,由于取的是矩形 ,使用以下的方式计算,如果取其他形状,就要修改这里
q1=p1;
q2=[p1(1) p1(2)+w];
q3=[p1(1)+h p1(2)];
q4=[p1(1)+h p1(2)+w];
%代入公式
B=[q1(1) q1(2) q2(1) q2(2) q3(1) q3(2) q4(1) q4(2)]';
A=[p1(1) p1(2) 1 0 0 0 -q1(1)*p1(1) -q1(1)*p1(2);
0 0 0 p1(1) p1(2) 1 -q1(2)*p1(1) -q1(2)*p1(2);
p2(1) p2(2) 1 0 0 0 -q2(1)*p2(1) -q2(1)*p2(2);
0 0 0 p2(1) p2(2) 1 -q2(2)*p2(1) -q2(2)*p2(2);
p3(1) p3(2) 1 0 0 0 -q3(1)*p3(1) -q3(1)*p3(2);
0 0 0 p3(1) p3(2) 1 -q3(2)*p3(1) -q3(2)*p3(2);
p4(1) p4(2) 1 0 0 0 -q4(1)*p4(1) -q4(1)*p4(2);
0 0 0 p4(1) p4(2) 1 -q4(2)*p4(1) -q4(2)*p4(2)];
%得到变换系数
tt=inv(A)*B;
%得到变换矩阵
perTra=[tt(1) tt(2) tt(3);
tt(4) tt(5) tt(6);
tt(7) tt(8) 1 ]
%首先要确定图像的边界,MATLAB中的数字是从1开始的,左上角的点应该是(1,1)而不是(0,0)
edge=perTra*[1 row 1 row;
1 1 col col;
1 1 1 1 ];
leftUp=[edge(1,1)/edge(3,1) edge(2,1)/edge(3,1)]
leftDown=[edge(1,2)/edge(3,2) edge(2,2)/edge(3,2)]
rightUp=[edge(1,3)/edge(3,3) edge(2,3)/edge(3,3)]
rightDown=[edge(1,4)/edge(3,4) edge(2,4)/edge(3,4)]
%变换后的图像的宽和高
height=round(max([leftUp(1) leftDown(1) rightUp(1) rightDown(1)])-min([leftUp(1) leftDown(1) rightUp(1) rightDown(1)]))
width=round(max([leftUp(2) leftDown(2) rightUp(2) rightDown(2)])-min([leftUp(2) leftDown(2) rightUp(2) rightDown(2)]))
imgOut=zeros(height,width);%定义新图像
upEdge=round(abs(min([leftUp(1) leftDown(1) rightUp(1) rightDown(1)])));
leftEdge=round(abs(min([leftUp(2) leftDown(2) rightUp(2) rightDown(2)])));
%变换后的图像与原图像可能会有很大的差异,从变换后的图像中,找所对应的原始图像的坐标的点
%这样做可以防止变换后的图像中出现空缺或多重的映射
invPerTra=inv(perTra);
for i = 1-upEdge:height-upEdge
for j = 1-leftEdge:width-leftEdge
originalPix=invPerTra*[i,j,1]';
originalPix(1)=originalPix(1)/originalPix(3);
originalPix(2)=originalPix(2)/originalPix(3);
if originalPix(1)>=1 && originalPix(2)>=1 && originalPix(1)<=row && originalPix(2)<=col
imgOut(i+upEdge,j+leftEdge)=img(round(originalPix(1)),round(originalPix(2))); %这里偷懒啦!使用最简单的最邻近插值
end
end
end
figure;
imshow(uint8(imgOut));
原图:
结果:
C语言版(为了保证移植性,没有使用函数库):
void PerspectiveTransformation(VERTEX LU,VERTEX RU,VERTEX LD,VERTEX RD)
{
int i,j,k,l,m;
int flag;
int wid,hig;
double width,high;
int dx,dy;
double P1[3][1]={1,1,1};
double PE[3][1]={0};
// double B[8][1];
// double A[8][8];
VERTEX LU1,LD1,RU1,RD1;//定义四个顶点
VERTEXINT LUn,LDn,RUn,RDn;//变换后的四个顶点
///乱给的数/
LU.x=46;
LU.y=24;
RU.x=470;
RU.y=48;
LD.x=15;
LD.y=278;
RD.x=470;
RD.y=318;
width=(sqrt((LU.x-RU.x)*(LU.x-RU.x)+(LU.y-RU.y)*(LU.y-RU.y))+0.5);//新矩形宽
high=(sqrt((LU.x-LD.x)*(LU.x-LD.x)+(LU.y-LD.y)*(LU.y-LD.y))+0.5);//新矩形高
printf("width:%f\n",width);
printf("high:%f\n",high);
变换后的顶点,方程右边的数/
LU1.x=LU.x;
LU1.y=LU.y;
RU1.x=LU.x+width;
RU1.y=LU.y;
LD1.x=LU.x;
LD1.y=LU.y+high;
RD1.x=LU.x+width;
RD1.y=LU.y+high;
{
double B[8][1]={LU1.x,LU1.y, RU1.x, RU1.y, LD1.x, LD1.y, RD1.x ,RD1.y};// %变换后的四个顶点,方程右边的值
double A[8][8]={
{LU.x , LU.y, 1 , 0 , 0 , 0, -LU1.x*LU.x, -LU1.x*LU.y},
{ 0 , 0 , 0 , LU.x, LU.y, 1, -LU1.y*LU.x, -LU1.y*LU.y},
{RU.x , RU.y, 1 , 0 , 0 , 0, -RU1.x*RU.x, -RU1.x*RU.y},
{ 0 , 0 , 0 , RU.x, RU.y, 1, -RU1.y*RU.x, -RU1.y*RU.y},
{LD.x , LD.y, 1 , 0 , 0 , 0, -LD1.x*LD.x, -LD1.x*LD.y},
{ 0 , 0 , 0 , LD.x, LD.y, 1, -LD1.y*LD.x, -LD1.y*LD.y},
{RD.x , RD.y, 1 , 0 , 0 , 0, -RD1.x*RD.x, -RD1.x*RD.y},
{ 0 , 0 , 0 , RD.x, RD.y, 1, -RD1.y*RD.x, -RD1.y*RD.y}};
/
inv_c(A);//求逆运算,可以参考我的另一篇boke
mulMatrix4B(invA,B,SIZE,SIZE,1,FA);
}
{
double pers[3][3]={{FA[3][0],FA[4][0],FA[5][0]},{FA[0][0],FA[1][0],FA[2][0]},{FA[6][0],FA[7][0],1}};
mulMatrix4P(pers,P1,3,3,1,PE); //复数运算
for(i=0;i<3;++i)
{
PE[i][0]/=(FA[6][0]+FA[7][0]+1);
//printf("LU:%0.31f\n\n",PE[i][0]);
}
LUn.x=(int)(PE[1][0]+0.5);
LUn.y=(int)(PE[0][0]+0.5);
printf("LUn.x:%d\nLUn.y:%d\n",LUn.x,LUn.y);
for(i=0;i<3;++i)
{
PE[i][0]=0;
}
P1[1][0]=IMGCOLSp;
mulMatrix4P(pers,P1,3,3,1,PE);
for(i=0;i<3;++i)
{
PE[i][0]/=(FA[6][0]+FA[7][0]*IMGCOLSp+1);
//printf("RU:%0.31f\n\n",PE[i][0]);
}
RUn.x=(int)(PE[1][0]+0.5);
RUn.y=(int)(PE[0][0]+0.5);
printf("RUn.x:%d\nRUn.y:%d\n",RUn.x,RUn.y);
for(i=0;i<3;++i)
{
PE[i][0]=0;
}
P1[0][0]=IMGROWSp;
P1[1][0]=1;
mulMatrix4P(pers,P1,3,3,1,PE);
for(i=0;i<3;++i)
{
PE[i][0]/=(FA[6][0]*IMGROWSp+FA[7][0]+1);
//printf("LD:%0.31f\n\n",PE[i][0]);
}
LDn.x=(int)(PE[1][0]+0.5);
LDn.y=(int)(PE[0][0]+0.5);
printf("LD.x:%d\nLD.y:%d\n",LDn.x,LDn.y);
for(i=0;i<3;++i)
{
PE[i][0]=0;
}
P1[1][0]=IMGCOLSp;
mulMatrix4P(pers,P1,3,3,1,PE);
for(i=0;i<3;++i)
{
PE[i][0]/=(FA[6][0]*IMGROWSp+FA[7][0]*IMGCOLSp+1);
}
RDn.x=(int)(PE[1][0]+0.5);
RDn.y=(int)(PE[0][0]+0.5);
printf("RD.x:%d\nRD.y:%d\n",RDn.x,RDn.y);
inv_cp(pers);
}
if(LUn.x>LDn.x)
{
k=LUn.x;
l=LDn.x;
}
else
{
k=LDn.x;
l=LUn.x;
}
if(k<RUn.x)
{
k=RUn.x;
}
if(k<RDn.x)
{
k=RDn.x;
}
if(l>RUn.x)
{
l=RUn.x;
}
if(l>RDn.x)
{
l=RDn.x;
}
//printf("big:%d\nsmall:%d\n",k,l);
wid=k-l;
dx=l>0?l:(-l);
//printf("wid:%d\n",wid);
if(LUn.y>LDn.y)
{
k=LUn.y;
l=LDn.y;
}
else
{
k=LDn.y;
l=LUn.y;
}
if(k<RUn.y)
{
k=RUn.y;
}
if(k<RDn.y)
{
k=RDn.y;
}
if(l>RUn.y)
{
l=RUn.y;
}
if(l>RDn.y)
{
l=RDn.y;
}
//printf("big:%d\nsmall:%d\n",k,l);
hig=k-l;
dy=l>0?l:(-l);
printf("hig:%d\n",hig);
printf("wid:%d\n",wid);
for(i=0;i<hig;++i)
{
for(j=0;j<wid;++j)
{
P1[0][0]=i;
P1[1][0]=j;
P1[2][0]=1;
mulMatrix4P(persOut,P1,3,3,1,PE);
{
double C[2][2]={FA[6][0]*PE[0][0]-1,FA[7][0]*PE[0][0],FA[6][0]*PE[1][0],FA[7][0]*PE[1][0]-1};
double CC[2][1]={-PE[0][0],-PE[1][0]};
inv_c2(C);
mulMatrix42(out2,CC,2,2,1,P2);
//printf("%f\n",C[1][1]);
}
// printf("%f\n",P2[0][0]);
if(P2[0][0]>=0.5 && P2[1][0]>=0.5 && P2[0][0]<=IMGROWSp && P2[1][0]<=IMGCOLSp)
{
imageOut[i][j]=BMP2.grayMatrix[(int)P2[0][0]][(int)P2[1][0]]; //最邻近插值
}
// {
// printf("%f\n",P2[0][0]);
// printf("%f\n",P2[1][0]);
// }
if(m<P2[0][0])
{
m=P2[0][0];
}
}
}
/
printf("max:%d\n",m);
bmpOutput8("yy.bmp",hig,wid);//保存为BMP文件
}
— 2022.11.17 —
理论部分这里直接贴图了哈。引用一篇论文的内容。
有同学对源码感兴趣,这里把文件直接贴出来了。多年前写的了,当时刚学编程,至于风格嘛,勿喷。
#include <stdio.h>
#include <math.h>
typedef struct {float re; float im;} COMPLEX;
typedef struct {double re; double im;} DCOMPLEX;
typedef struct {double x;double y;} VERTEX;
typedef struct {int x;int y;} VERTEXINT;
#define PI 3.14159265358979643864246
#define PI2 6.28318530717958647728492 /* 2.0 * PI */
#define FFT_FOR 0//正变换标志
#define FFT_INV 1//逆变换标志
#define IMGROWS 512 //行数
#define IMGCOLS 512 //列数
#define IMGROWSp 512 //p行数
#define IMGCOLSp 512 //p列数
#define BLOCK 64
#define SIZE 8//透视变换用到的方阵
struct BMPHEADER{//小端方式存储,麻烦!!
unsigned char ucBuf[54];//文件头信息存放的数组
unsigned int intBuf[54];
unsigned int bmpSize;//文件大小
unsigned int bmpOffset;//数据偏移量
unsigned char infoSize;//位图信息头的大小
int cols;//列数 图像的宽度
int rows;//行数 图像的高度
unsigned int pointSize;//位深度 每个像素的大小
char compress; //图像是否压缩
unsigned char RGBquad[1024];//调色板
unsigned char mateix[IMGROWS][IMGCOLS*3];//存放图片数据的数组
unsigned char grayMatrix[IMGROWS][IMGCOLS];
}BMPHEAD, BMP1, BMP2,OUT;
unsigned char imageShow[IMGROWS][IMGCOLS]={0};
int deviationX=0;
int deviationY=0;
int startX=0;
int startY=0;
int disX=0;
int disY=0;
int deviationLUX=0;
int deviationLUY=0;
int deviationLDX=0;
int deviationLDY=0;
int deviationRUX=0;
int deviationRUY=0;
int deviationRDX=0;
int deviationRDY=0;
unsigned char imageOut[IMGROWS*2][IMGCOLS*2]={0};//定义时取极限,写入图片时再作调整
unsigned char imageOut24[IMGROWS*2][IMGCOLS*6]={0};
///
int bmpOutput8(FILE * fpTarget,unsigned int row,unsigned int col)//只用于平移的图像
{
int i=0,j=0;
OUT=BMP1;
/// 必须对齐,不然会跪
OUT.rows=(row+3)/4*4;
OUT.cols=(col+3)/4*4;
fpTarget = fopen (fpTarget, "wb") ;
OUT.ucBuf[0]=0x42;
OUT.ucBuf[1]=0x4d;
OUT.bmpSize=OUT.rows*OUT.cols+1024+54;
OUT.ucBuf[2]=(unsigned char)OUT.bmpSize;
OUT.ucBuf[3]=(unsigned char)(OUT.bmpSize>>8);
OUT.ucBuf[4]=(unsigned char)(OUT.bmpSize>>16);
OUT.ucBuf[5]=(unsigned char)(OUT.bmpSize>>24);
OUT.ucBuf[6]=0;
OUT.ucBuf[7]=0;
OUT.ucBuf[8]=0;
OUT.ucBuf[9]=0;
OUT.bmpOffset=1078;
OUT.ucBuf[10]=(unsigned char)OUT.bmpOffset;
OUT.ucBuf[11]=(unsigned char)(OUT.bmpOffset>>8);
OUT.ucBuf[12]=(unsigned char)(OUT.bmpOffset>>16);
OUT.ucBuf[13]=(unsigned char)(OUT.bmpOffset>>24);
OUT.infoSize=1064;
OUT.ucBuf[14]=(unsigned char)OUT.infoSize;
OUT.ucBuf[15]=(unsigned char)(OUT.infoSize>>8);
OUT.ucBuf[16]=(unsigned char)(OUT.infoSize>>16);
OUT.ucBuf[17]=(unsigned char)(OUT.infoSize>>24);
OUT.ucBuf[18]=(unsigned char)OUT.cols;
OUT.ucBuf[19]=(unsigned char)(OUT.cols>>8);
OUT.ucBuf[20]=(unsigned char)(OUT.cols>>16);
OUT.ucBuf[21]=(unsigned char)(OUT.cols>>24);
OUT.ucBuf[22]=(unsigned char)OUT.rows;
OUT.ucBuf[23]=(unsigned char)(OUT.rows>>8);
OUT.ucBuf[24]=(unsigned char)(OUT.rows>>16);
OUT.ucBuf[25]=(unsigned char)(OUT.rows>>24);
OUT.ucBuf[26]=0x01;
OUT.ucBuf[27]=0;
OUT.ucBuf[28]=0x08;
OUT.ucBuf[29]=0x00;
OUT.ucBuf[30]=0;
OUT.ucBuf[31]=0;
OUT.ucBuf[32]=0;
OUT.ucBuf[33]=0;
OUT.ucBuf[34]=0;
OUT.ucBuf[35]=0;
OUT.ucBuf[36]=0;
//OUT.cols=col;
///
OUT.ucBuf[37]=0;
OUT.ucBuf[38]=0;
OUT.ucBuf[39]=0;
OUT.ucBuf[40]=0;
OUT.ucBuf[41]=0;
OUT.ucBuf[42]=0;
OUT.ucBuf[43]=0;
OUT.ucBuf[44]=0;
OUT.ucBuf[45]=0;
fwrite (OUT.ucBuf, 1, sizeof(OUT.ucBuf), fpTarget);
for(i=0;i<256;i++)//生成调色板
{
OUT.RGBquad[4*i]=i;
OUT.RGBquad[4*i+1]=i;
OUT.RGBquad[4*i+2]=i;
OUT.RGBquad[4*i+3]=0;
}
fwrite (OUT.RGBquad, 1, sizeof(OUT.RGBquad), fpTarget);
for(i=OUT.rows-1;i>=0;--i)
{
// for(j=0;j<OUT.cols;++j)
// {
//printf("i:%d j:%d\n",i,j);
fwrite (imageOut[i], 1 ,OUT.cols, fpTarget);
//}
}
fclose (fpTarget) ;
return 1;
}
int bmp2matrix(FILE * fpPhoto)
{
int i=0,j=0;
int *point;
//int **arr;
// int zz[32000];
unsigned int readResult;
fread(&BMPHEAD.ucBuf,1,54,fpPhoto);
for(i=0;i<54;++i)
{
printf("%d %X\n",i,BMPHEAD.ucBuf[i]);
}
if(BMPHEAD.ucBuf[0]==0x42&&BMPHEAD.ucBuf[1]==0x4d)
{
for(i=2;i<54;i++)//数据类型与存储空间
{
BMPHEAD.intBuf[i]=BMPHEAD.ucBuf[i];
// printf("%d\n",BMPHEAD.intBuf[i]);
}
BMPHEAD.bmpSize=(BMPHEAD.intBuf[5]<<24)+(BMPHEAD.intBuf[4]<<16)+(BMPHEAD.intBuf[3]<<8)+(BMPHEAD.intBuf[2]);//计算BMP图片的大小,左移符号优先级低!!
//华丽的分隔符,主角登场!/
BMPHEAD.bmpOffset=(BMPHEAD.intBuf[13]<<24)+(BMPHEAD.intBuf[12]<<16)+(BMPHEAD.intBuf[11]<<8)+(BMPHEAD.intBuf[10]);//a-d 获得数据信息的偏移量
// printf("offset: %d \n",BMPHEAD.bmpOffset);
///
BMPHEAD.infoSize=(BMPHEAD.intBuf[17]<<24)+(BMPHEAD.intBuf[16]<<16)+(BMPHEAD.intBuf[15]<<8)+(BMPHEAD.intBuf[14]);//0e-11 位图信息头的大小
// printf("infoSize: %d \n",BMPHEAD.infoSize);
BMPHEAD.cols=(BMPHEAD.intBuf[21]<<24)+(BMPHEAD.intBuf[20]<<16)+(BMPHEAD.intBuf[19]<<8)+(BMPHEAD.intBuf[18]);//12-15 位图的宽度
printf("cols: %d \n",BMPHEAD.cols);
BMPHEAD.rows=(BMPHEAD.intBuf[25]<<24)+(BMPHEAD.intBuf[24]<<16)+(BMPHEAD.intBuf[23]<<8)+(BMPHEAD.intBuf[22]);//16-19 位图的高度
// printf("rows: %d \n",BMPHEAD.rows);
BMPHEAD.pointSize=(BMPHEAD.intBuf[29]<<8)+(BMPHEAD.intBuf[28]);//1c-1d 位深度
// printf("pointSize: %d \n",BMPHEAD.pointSize);
if(BMPHEAD.pointSize==24)
{
//这是24位图,没有调色板
}
switch (BMPHEAD.pointSize)
{
case 1:
printf("位深度出错!!!!\n");
//二值化 8个像素占1个字节
break;
case 4:
printf("位深度出错!!!!\n");
// 2个像素占1个字节
break;
case 8:
printf("位深度出错!!!!\n");
//灰度图 1个像素占1个字节
break;
case 24:
printf("BMP__24\n");
//24位位图 1个像素占3个字节 没有调色板,暂时只处理这种图,其他的以后优化时再说
break;
case 32:
printf("位深度出错!!!!\n");
//32位位图 没有调色板,多了一个Alpha通道
break;
default:printf("位深度出错!!!!\n");break;
}
if(0==(BMPHEAD.compress=(BMPHEAD.intBuf[33]<<24)+(BMPHEAD.intBuf[32]<<16)+(BMPHEAD.intBuf[31]<<8)+(BMPHEAD.intBuf[30])))//1e-21 为0表示不压缩
{
printf("no compress\n");
}
// point=(int *)malloc(BMPHEAD.bmpSize);
// if(point==NULL)
// {
// printf("内存分配失败!!!!!!");
// }
fread(&point,1,BMPHEAD.bmpOffset-54,fpPhoto);
for(i=BMPHEAD.rows-1;i>=0;i--)
{
readResult=fread(&BMPHEAD.mateix[i],1,3*BMPHEAD.cols,fpPhoto);//fread读过之后文件指针移动,就是指读过的东西,就读没了!!!!!hehehehehehehehehe
//24位图,三个通道,列数X3
}//将图片读到数组mateix中,在数据区左下是第一个点 byte:B byte:G byte:R
//printf("成功!!!!!!");
fclose (fpPhoto) ;
return 1;
}
else
{
printf("这不是BMP!!");
fclose (fpPhoto) ;
return 0; //不正常 返回0
}
}
DCOMPLEX first[IMGROWS][IMGCOLS];
DCOMPLEX second[IMGROWS][IMGCOLS];
DCOMPLEX relevant[IMGROWS][IMGCOLS];
DCOMPLEX overlap1[IMGROWS][IMGCOLS]={0};
DCOMPLEX overlap2[IMGROWS][IMGCOLS]={0};
DCOMPLEX overlap1LU[BLOCK][BLOCK]={0};
DCOMPLEX overlap1LD[BLOCK][BLOCK]={0};
DCOMPLEX overlap1RU[BLOCK][BLOCK]={0};
DCOMPLEX overlap1RD[BLOCK][BLOCK]={0};
DCOMPLEX overlap2LU[BLOCK][BLOCK]={0};
DCOMPLEX overlap2LD[BLOCK][BLOCK]={0};
DCOMPLEX overlap2RU[BLOCK][BLOCK]={0};
DCOMPLEX overlap2RD[BLOCK][BLOCK]={0};
DCOMPLEX overlaprele[BLOCK][BLOCK]={0};
/*/
函数名:void mulMatrix(double matrix1[SIZE][SIZE],double matrix2[SIZE][SIZE],int high1,int weight,int weight2,double mulMatrixOut[SIZE][SIZE])
输入:依次是 左矩阵,右矩阵,左矩阵高度,左矩阵宽度,右矩阵宽度,输出矩阵
输出:
功能:矩阵乘法
作者:dogdng
*/
void mulMatrix(double matrix1[SIZE][SIZE],double matrix2[SIZE][SIZE],int high1,int weight,int weight2,double mulMatrixOut[SIZE][SIZE])
{
int i,j,k;
for(i=0;i<high1;++i)
{
for(j=0;j<weight2;j++)
{
for(k=0;k<weight;++k)
{
mulMatrixOut[i][j]+=matrix1[i][k]*matrix2[k][j];
}
}
}
}
void mulMatrix3(double matrix1[3][3],double matrix2[3][3],int high1,int weight,int weight2,double mulMatrixOut[3][3])
{
int i,j,k;
for(i=0;i<high1;++i)
{
for(j=0;j<weight2;j++)
{
for(k=0;k<weight;++k)
{
mulMatrixOut[i][j]+=matrix1[i][k]*matrix2[k][j];
}
}
}
}
void mulMatrix4B(double matrix1[SIZE][SIZE],double matrix2[SIZE][1],int high1,int weight,int weight2,double mulMatrixOut[SIZE][1])
{
int i,j,k;
for(i=0;i<SIZE;i++)
for(j=0;j<SIZE;j++)
mulMatrixOut[i][j]=0;
for(i=0;i<high1;++i)
{
for(j=0;j<weight2;j++)
{
for(k=0;k<weight;++k)
{
mulMatrixOut[i][j]+=matrix1[i][k]*matrix2[k][j];
}
}
//printf("i am ok\n");
}
}
void mulMatrix4P(double matrix1[3][3],double matrix2[3][1],int high1,int weight,int weight2,double mulMatrixOut[3][1])
{
int i,j,k;
for(i=0;i<3;++i)
mulMatrixOut[i][0]=0;
for(i=0;i<high1;++i)
{
for(j=0;j<weight2;j++)
{
for(k=0;k<weight;++k)
{
mulMatrixOut[i][j]+=matrix1[i][k]*matrix2[k][j];
}
}
//printf("i am ok\n");
}
}
void mulMatrix42(double matrix1[2][2],double matrix2[2][1],int high1,int weight,int weight2,double mulMatrixOut[2][1])
{
mulMatrixOut[0][0]=matrix1[0][0]*matrix2[0][0]+matrix1[0][1]*matrix2[1][0];
mulMatrixOut[1][0]=matrix1[1][0]*matrix2[0][0]+matrix1[1][1]*matrix2[1][0];
}
/*/
函数名:void inv_c(double A[SIZE][SIZE])
输入:
输出:
功能:求矩阵的逆 pure C language
首先对矩阵进行QR分解之后求上三角矩阵R的逆阵最后A-1=QH*R-1,得到A的逆阵。
作者:dogdng
*/
double invA[SIZE][SIZE]={0};//A的逆矩阵,最终的结果
double FA[SIZE][1]={0};
void inv_c(double A[SIZE][SIZE])
{
int i;//数组 行
int j;//数组 列
int k;//代表B的角标
int l;//数组 列
double dev;
double numb;//计算的中间变量
double numerator,denominator;
double ratio;
double b[SIZE][SIZE]={0};//应该读作“贝尔塔”,注释中用B表示
double t[SIZE][SIZE]={0};//求和的那项
double Q[SIZE][SIZE]={0};//正交矩阵
double QH[SIZE][SIZE]={0};//正交矩阵的转置共轭
double R[SIZE][SIZE]={0};//
double invR[SIZE][SIZE]={0};//R的逆矩阵
//={0};//
double matrixR1[SIZE][SIZE]={0};
double matrixR2[SIZE][SIZE]={0};
/*求B*/
for(i=0;i<SIZE;++i)
{
for(j=0;j<SIZE;++j)
{
b[j][i]=A[j][i];
}
for(k=0;k<i;++k)
{
if(i)
{
numerator=0.0;
denominator=0.0;
for(l=0;l<SIZE;++l)
{
numerator+=A[l][i]*b[l][k];
denominator+=b[l][k]*b[l][k];
}
dev=numerator/denominator;
t[k][i]=dev;
for(j=0;j<SIZE;++j)
{
b[j][i]-=t[k][i]*b[j][k];//t init =0 !!!
}
}
}
}
///对B单位化,得到正交矩阵Q矩阵
for(i=0;i<SIZE;++i)
{
numb=0.0;
for(j=0;j<SIZE;++j)
{
numb+=(b[j][i]*b[j][i]);
}
dev=sqrt(numb);
for(j=0;j<SIZE;++j)
{
Q[j][i]=b[j][i]/dev;
}
matrixR1[i][i]=dev;
}
//求上三角R阵///
for(i=0;i<SIZE;++i)
{
for(j=0;j<SIZE;++j)
{
if(j<i)
{
matrixR2[j][i]=t[j][i];
}
else if(j==i)
{
matrixR2[j][i]=1;
}
else
{
matrixR2[j][i]=0;
}
}
}
mulMatrix(matrixR1,matrixR2,SIZE,SIZE,SIZE,R);
///QR分解完毕//
printf("QR分解:\n");
printf("Q=\n");
for(i=0;i<SIZE;++i)
{
for(j=0;j<SIZE;++j)
{
printf("%2.4f ",Q[i][j]);
//
}
printf("\n");
}
printf("R=\n");
for(i=0;i<SIZE;++i)
{
for(j=0;j<SIZE;++j)
{
printf("%2.4f ",R[i][j]);
//
}
printf("\n");
}
//求R的逆阵//
for(i=SIZE-1;i>=0;--i)
{
invR[i][i]=1/R[i][i];
//R[i][i]=invR[i][i];
if(i!=(SIZE-1))//向右
{
for(j=i+1;j<SIZE;++j)
{
invR[i][j]=invR[i][j]*invR[i][i];
R[i][j]=R[i][j]*invR[i][i];
}
}
if(i)//向上
{
for(j=i-1;j>=0;--j)
{
ratio=R[j][i];
for(k=i;k<SIZE;++k)
{
invR[j][k]-=ratio*invR[i][k];
R[j][k]-=ratio*R[i][k];
}
}
}
}
///
printf("inv(R)=\n");
for(i=0;i<SIZE;++i)
{
for(j=0;j<SIZE;++j)
{
printf(" %2.4f ",invR[i][j]);
//
}
printf("\n");
}
//结果和MATLAB差一个负号,神马鬼????????/
//求QH//
for(i=0;i<SIZE;++i)//实矩阵就是转置
{
for(j=0;j<SIZE;++j)
{
QH[i][j]=Q[j][i];
}
}
///求A的逆阵invA/
mulMatrix(invR,QH,SIZE,SIZE,SIZE,invA);
printf("inv(A)=\n");
for(i=0;i<SIZE;++i)
{
for(j=0;j<SIZE;++j)
{
printf(" %2.4f ",invA[i][j]);
//
}
printf("\n");
}
///结果与MATLAB的结果在千分位后有出入,但是负号都是对的^v^///
}
double persOut[3][3]={0};
void inv_cp(double A[3][3])
{
int i;//数组 行
int j;//数组 列
int k;//代表B的角标
int l;//数组 列
double dev;
double numb;//计算的中间变量
double numerator,denominator;
double ratio;
double b[3][3]={0};//应该读作“贝尔塔”,注释中用B表示
double t[3][3]={0};//求和的那项
double Q[3][3]={0};//正交矩阵
double QH[3][3]={0};//正交矩阵的转置共轭
double R[3][3]={0};//
double invR[3][3]={0};//R的逆矩阵
//={0};//
double matrixR1[3][3]={0};
double matrixR2[3][3]={0};
/*求B*/
for(i=0;i<3;++i)
{
for(j=0;j<3;++j)
{
b[j][i]=A[j][i];
}
for(k=0;k<i;++k)
{
if(i)
{
numerator=0.0;
denominator=0.0;
for(l=0;l<3;++l)
{
numerator+=A[l][i]*b[l][k];
denominator+=b[l][k]*b[l][k];
}
dev=numerator/denominator;
t[k][i]=dev;
for(j=0;j<3;++j)
{
b[j][i]-=t[k][i]*b[j][k];//t init =0 !!!
}
}
}
}
///对B单位化,得到正交矩阵Q矩阵
for(i=0;i<3;++i)
{
numb=0.0;
for(j=0;j<3;++j)
{
numb+=(b[j][i]*b[j][i]);
}
dev=sqrt(numb);
for(j=0;j<3;++j)
{
Q[j][i]=b[j][i]/dev;
}
matrixR1[i][i]=dev;
}
//求上三角R阵///
for(i=0;i<3;++i)
{
for(j=0;j<3;++j)
{
if(j<i)
{
matrixR2[j][i]=t[j][i];
}
else if(j==i)
{
matrixR2[j][i]=1;
}
else
{
matrixR2[j][i]=0;
}
}
}
mulMatrix3(matrixR1,matrixR2,3,3,3,R);
///QR分解完毕//
printf("QR分解:\n");
printf("Q=\n");
for(i=0;i<3;++i)
{
for(j=0;j<3;++j)
{
printf("%2.4f ",Q[i][j]);
//
}
printf("\n");
}
printf("R=\n");
for(i=0;i<3;++i)
{
for(j=0;j<3;++j)
{
printf("%2.4f ",R[i][j]);
//
}
printf("\n");
}
//求R的逆阵//
for(i=3-1;i>=0;--i)
{
invR[i][i]=1/R[i][i];
//R[i][i]=invR[i][i];
if(i!=(3-1))//向右
{
for(j=i+1;j<3;++j)
{
invR[i][j]=invR[i][j]*invR[i][i];
R[i][j]=R[i][j]*invR[i][i];
}
}
if(i)//向上
{
for(j=i-1;j>=0;--j)
{
ratio=R[j][i];
for(k=i;k<3;++k)
{
invR[j][k]-=ratio*invR[i][k];
R[j][k]-=ratio*R[i][k];
}
}
}
}
///
printf("inv(R)=\n");
for(i=0;i<3;++i)
{
for(j=0;j<3;++j)
{
printf(" %2.4f ",invR[i][j]);
//
}
printf("\n");
}
//结果和MATLAB差一个负号,神马鬼????????/
//求QH//
for(i=0;i<3;++i)//实矩阵就是转置
{
for(j=0;j<3;++j)
{
QH[i][j]=Q[j][i];
}
}
///求A的逆阵invA/
mulMatrix3(invR,QH,3,3,3,persOut);
}
double out2[2][2]={0};
void inv_c2(double A[2][2])
{
double i;
i=A[0][0]*A[1][1]-A[0][1]*A[1][0];
out2[0][0]=A[1][1]/i;
out2[0][1]=-A[0][1]/i;
out2[1][0]=-A[1][0]/i;
out2[1][1]=A[0][0]/i;
}
double P2[2][1]={0};
int main()
{
int i,j,k,l,m,n;
int flag;
int wid,hig;
double sm;//计算revelant的中间变量
DCOMPLEX max,min;//归一化使用的极值
double maxmax=0,minmin=0;
double modMax=0;
double width,high;
int dx,dy;
double P1[3][1]={1,1,1};
double PE[3][1]={0};
// double B[8][1];
// double A[8][8];
VERTEX LU,LD,RU,RD;//定义四个顶点
VERTEX LU1,LD1,RU1,RD1;//定义四个顶点
VERTEXINT LUn,LDn,RUn,RDn;//变换后的四个顶点
FILE * fpPhoto1,* fpPhoto2;
// fpPhoto1 = fopen ("ru.bmp", "rb") ;//输入图片1为512*512
// fpPhoto2 = fopen ("ld.bmp", "rb") ;//输入图片2为512*512
fpPhoto1 = fopen ("11.bmp", "rb") ;//输入图片1为512*512
fpPhoto2 = fopen ("12.bmp", "rb") ;//输入图片2为512*512
if (!fpPhoto1||!fpPhoto2)//文件打开失败则退出
{
printf ("open file faild!\n");
return -1 ;
}/*if (!fpPhoto)*/
if(!bmp2matrix(fpPhoto1))
{
printf ("open BMP1 faild!\n");
return -1 ;
}
BMP1=BMPHEAD;
if(!bmp2matrix(fpPhoto2))
{
printf ("open BMP2 faild!\n");
return -1 ;
}
BMP2=BMPHEAD;
//彩色信息转为灰度存到复数中//
for(i=0;i<IMGROWS;i++)
{
for(j=0;j<IMGCOLS;j++)
{
BMP1.grayMatrix[i][j]=(BMP1.mateix[i][j*3]*11+ BMP1.mateix[i][j*3+1]*59+BMP1.mateix[i][j*3+2] *30)/100;
BMP2.grayMatrix[i][j]=(BMP2.mateix[i][j*3]*11+ BMP2.mateix[i][j*3+1]*59+BMP2.mateix[i][j*3+2] *30)/100;
}
}
///
{
///乱给的数/
LU.x=61;
LU.y=73;
RU.x=60;
RU.y=357;
LD.x=337;
LD.y=127;
RD.x=330;
RD.y=490;
width=(sqrt((LU.x-RU.x)*(LU.x-RU.x)+(LU.y-RU.y)*(LU.y-RU.y))+0.5);//新矩形宽
high=(sqrt((LU.x-LD.x)*(LU.x-LD.x)+(LU.y-LD.y)*(LU.y-LD.y))+0.5);//新矩形高
printf("width:%f\n",width);
printf("high:%f\n",high);
//变换后的顶点,方程右边的数/
LU1.x=LU.x;
LU1.y=LU.y;
RU1.x=LU.x+width;
RU1.y=LU.y;
LD1.x=LU.x;
LD1.y=LU.y+high;
RD1.x=LU.x+width;
RD1.y=LU.y+high;
{
double B[8][1]={LU1.x,LU1.y, RU1.x, RU1.y, LD1.x, LD1.y, RD1.x ,RD1.y};// %变换后的四个顶点,方程右边的值
double A[8][8]={
{LU.x , LU.y, 1 , 0 , 0 , 0, -LU1.x*LU.x, -LU1.x*LU.y},
{ 0 , 0 , 0 , LU.x, LU.y, 1, -LU1.y*LU.x, -LU1.y*LU.y},
{RU.x , RU.y, 1 , 0 , 0 , 0, -RU1.x*RU.x, -RU1.x*RU.y},
{ 0 , 0 , 0 , RU.x, RU.y, 1, -RU1.y*RU.x, -RU1.y*RU.y},
{LD.x , LD.y, 1 , 0 , 0 , 0, -LD1.x*LD.x, -LD1.x*LD.y},
{ 0 , 0 , 0 , LD.x, LD.y, 1, -LD1.y*LD.x, -LD1.y*LD.y},
{RD.x , RD.y, 1 , 0 , 0 , 0, -RD1.x*RD.x, -RD1.x*RD.y},
{ 0 , 0 , 0 , RD.x, RD.y, 1, -RD1.y*RD.x, -RD1.y*RD.y}};
//
printf("A=\n");
for(i=0;i<SIZE;++i)
{
for(j=0;j<SIZE;++j)
{
printf(" %2.4f ",A[i][j]);
//
}
printf("\n");
}
printf("\n");
printf("\n");
printf("\n");
for(i=0;i<SIZE;++i)
{
printf("%0.31f\n\n",B[i][0]);
}
printf("\n");
printf("\n");
printf("\n");
inv_c(A);
mulMatrix4B(invA,B,SIZE,SIZE,1,FA);
}
}
//FA is right!!!!!!!!!
{
double pers[3][3]={{FA[3][0],FA[4][0],FA[5][0]},{FA[0][0],FA[1][0],FA[2][0]},{FA[6][0],FA[7][0],1}};
mulMatrix4P(pers,P1,3,3,1,PE);
for(i=0;i<3;++i)
{
PE[i][0]/=(FA[6][0]+FA[7][0]+1);
//printf("LU:%0.31f\n\n",PE[i][0]);
}
LUn.x=(int)(PE[1][0]+0.5);
LUn.y=(int)(PE[0][0]+0.5);
printf("LUn.x:%d\nLUn.y:%d\n",LUn.x,LUn.y);
for(i=0;i<3;++i)
{
PE[i][0]=0;
}
P1[1][0]=IMGCOLSp;
mulMatrix4P(pers,P1,3,3,1,PE);
for(i=0;i<3;++i)
{
PE[i][0]/=(FA[6][0]+FA[7][0]*IMGCOLSp+1);
//printf("RU:%0.31f\n\n",PE[i][0]);
}
RUn.x=(int)(PE[1][0]+0.5);
RUn.y=(int)(PE[0][0]+0.5);
printf("RUn.x:%d\nRUn.y:%d\n",RUn.x,RUn.y);
for(i=0;i<3;++i)
{
PE[i][0]=0;
}
P1[0][0]=IMGROWSp;
P1[1][0]=1;
mulMatrix4P(pers,P1,3,3,1,PE);
for(i=0;i<3;++i)
{
PE[i][0]/=(FA[6][0]*IMGROWSp+FA[7][0]+1);
//printf("LD:%0.31f\n\n",PE[i][0]);
}
LDn.x=(int)(PE[1][0]+0.5);
LDn.y=(int)(PE[0][0]+0.5);
printf("LD.x:%d\nLD.y:%d\n",LDn.x,LDn.y);
for(i=0;i<3;++i)
{
PE[i][0]=0;
}
P1[1][0]=IMGCOLSp;
mulMatrix4P(pers,P1,3,3,1,PE);
for(i=0;i<3;++i)
{
PE[i][0]/=(FA[6][0]*IMGROWSp+FA[7][0]*IMGCOLSp+1);
}
RDn.x=(int)(PE[1][0]+0.5);
RDn.y=(int)(PE[0][0]+0.5);
printf("RD.x:%d\nRD.y:%d\n",RDn.x,RDn.y);
inv_cp(pers);
}
if(LUn.x>LDn.x)
{
k=LUn.x;
l=LDn.x;
}
else
{
k=LDn.x;
l=LUn.x;
}
if(k<RUn.x)
{
k=RUn.x;
}
if(k<RDn.x)
{
k=RDn.x;
}
if(l>RUn.x)
{
l=RUn.x;
}
if(l>RDn.x)
{
l=RDn.x;
}
//printf("big:%d\nsmall:%d\n",k,l);
wid=k-l;
dx=l>0?l:(-l);
//printf("wid:%d\n",wid);
if(LUn.y>LDn.y)
{
k=LUn.y;
l=LDn.y;
}
else
{
k=LDn.y;
l=LUn.y;
}
if(k<RUn.y)
{
k=RUn.y;
}
if(k<RDn.y)
{
k=RDn.y;
}
if(l>RUn.y)
{
l=RUn.y;
}
if(l>RDn.y)
{
l=RDn.y;
}
//printf("big:%d\nsmall:%d\n",k,l);
hig=k-l;
dy=l>0?l:(-l);
printf("hig:%d\n",hig);
printf("wid:%d\n",wid);
for(i=0;i<hig;++i)
{
for(j=0;j<wid;++j)
{
P1[0][0]=i;
P1[1][0]=j;
P1[2][0]=1;
mulMatrix4P(persOut,P1,3,3,1,PE);
{
double C[2][2]={FA[6][0]*PE[0][0]-1,FA[7][0]*PE[0][0],FA[6][0]*PE[1][0],FA[7][0]*PE[1][0]-1};
double CC[2][1]={-PE[0][0],-PE[1][0]};
inv_c2(C);
mulMatrix42(out2,CC,2,2,1,P2);
//printf("%f\n",C[1][1]);
}
// printf("%f\n",P2[0][0]);
if(P2[0][0]>=0.5 && P2[1][0]>=0.5 && P2[0][0]<=IMGROWSp && P2[1][0]<=IMGCOLSp)
{
imageOut[i][j]=BMP2.grayMatrix[(int)P2[0][0]][(int)P2[1][0]]; //最邻近插值,也可以用双线性或双立方插值
}
// {
// printf("%f\n",P2[0][0]);
// printf("%f\n",P2[1][0]);
// }
if(m<P2[0][0])
{
m=P2[0][0];
}
}
}
printf("max:%d\n",m);
bmpOutput8("yy.bmp",hig,wid);
printf("Vectory!!\n");
fclose (fpPhoto1);
fclose (fpPhoto2);
return 0;
}